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ABSTRACT 

We present a comparison of 14 galaxy formation models: 12 different semi-analytical 
models and 2 halo-occupation distribution models for galaxy formation based upon the same 
cosmological simulation and merger tree information derived from it. 

The participating codes have proven to be very successful in their own right but they 
have all been calibrated independently using various observational data sets, stellar models, 
and merger trees. In this paper we apply them without recalibration and this leads to a wide 
variety of predictions for the stellar mass function, specific star formation rates, stellar-to- 
halo mass ratios, and the abundance of orphan galaxies. The scatter is much larger than seen 
in previous comparison studies primarily because the codes have been used outside of their 
native environment within which they are well tested and calibrated. 

The purpose of the ‘nIFTy comparison of galaxy formation models’ is to bring together as 
many different galaxy formation modellers as possible and to investigate a common approach 
to model calibration. This paper provides a unified description for all participating models 
and presents the initial, uncalibrated comparison as a baseline for our future studies where we 
will develop a common calibration framework and address the extent to which that reduces 
the scatter in the model predictions seen here. 

Key words: methods: iV-body simulations - galaxies: haloes - galaxies: evolution - cosmol¬ 
ogy: theory - dark matter 


1 INTRODUCTION 

Understanding the formation and evolution of galaxies within a 
self-consistent cosmological context is one of the outstanding and 
most challenging topics of astrophysics and cosmology. Over the 
last few decades great strides forward have been made along two 
distinct lines: on the one hand, through directly accounting for the 
baryonic component (gas, stars, supermassive black holes, etc.) in 
cosmological simulations that include hydrodynamics and gravity 
and on the other hand, through a procedure known as semi-analytic 
modelling (SAM), in which a statistical estimate of the distribution 
of dark matter haloes and their merger history - either coming from 
cosmological simulations or extended Press-Schechter/Lagrangian 
methods - is combined with simplified yet physically motivated 
prescriptions to estimate the distribution of the physical proper¬ 
ties of galaxies. To date the vast computational challenge related 
to simulating the baryonic component has made it impractical for a 
general adoption of the former approach within the large volumes 
necessary for galaxy surveys and hence a lot of effort has been de¬ 
voted to the latter SAM strategy. Further, most of the modelling of 
sub-grid physics in hydrodynamical simulations relies on schemes 
akin to the ones used in the SAM approach, and hence it is more 
effective to apply them in post-processing to a simulation where pa¬ 
rameter scans are then less costly. Further, over this period we have 
not only witnessed significant advances in simulation techniques 
(e.g. Vogelsberger et al. 2014; Schaye et al. 2015), but the original 
ideas for SAMs (cf. White & Rees 1978) have undergone substan¬ 
tial refinement too (e.g. Cole 1991; White & Frenk 1991; Lacey 
& Silk 1991; Blanchard et al. 1992; Kauffmann et al. 1993, and in 
particular all the people and methods detailed below in Section 3). 

Some SAMs simply rely on analytical forms for the under¬ 
lying merger trees based upon (conditional) mass functions from 
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Press & Schechter (1974) or Extended Press-Schechter (Bond et al. 
1991) as, for instance, described in Somerville & Kolatt (1999). 
Other codes take as input halo merger trees derived from cosmolog¬ 
ical W-body simulations (see Lacey & Cole 1993; Roukema et al. 
1997, for the historical origin of both techniques). While the former 
remain a critical and powerful approach, advances in computing 
power (especially for dark matter only simulations) have shifted the 
focus of SAM developers towards the utilization of W-body merger 
trees as input to their models as they more reliably capture non¬ 
linear structure formation. For a recent comparison of merger tree 
construction methods, the influence of the underlying halo finder 
and the impact for (a particular) SAM we refer the reader to re¬ 
sults coming out of a previous comparison project ‘Sussing Merger 
Trees’^ (Srisawat et al. 2013; Avila et al. 2014; Lee et al. 2014): 
Lee et al. (2014), for instance, have shown that SAM parameters 
can be re-tuned to overcome differences between different merger 
trees; something of great relevance for the work presented here, as 
we will see later. 

As well as considering the influence of the halo finder and 
merger tree construction on any galaxy formation model, it is also 
important to consider the different semi-analytical techniques and 
methods themselves. Where this strategy has been used it has thus 
far focused primarily on comparing the physical details of the 
model. For instance, Somerville & Primack (1999) (as well as Lu 
et al. 2011) implemented various physical prescriptions into a sin¬ 
gle code; by design this tested the underlying physical assumptions 
and principles rather than for any code-to-code (dis-)similarities. 
Fontanot et al. (2009) and Kimm et al. (2009) compared differ¬ 
ent SAMs, but without using the same merger trees as an input. 
Similarly, Contreras et al. (2013) compared Durham and Munich 
SAMs for the Millennium simulation (Springel et al. 2005) mea¬ 
suring and comparing the halo occupation distribution found within 
them. Different channels for bulge formation in two distinct SAMs 
have been compared by De Lucia et al. (2011) and Fontanot et al. 
(2011). Fontanot et al. (2012) compared the predictions of three 

^ http :/ /popia . ft . uam.es/SussingMergerTrees 
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different SAMs (the same ones as used already in Fontanot et al. 
2009; Kimm et al. 2009) for the star formation rate function to ob¬ 
servations. 

The first move towards using identical inputs was undertaken 
by Maccio et al. (2010) where three SAMs were compared, this 
time using the same merger tree for all of them. However, the em¬ 
phasis was on studying (four) Milky Way-sized dark matter haloes. 
Dfaz-Gimenez & Mamon (2010) have analyzed the properties of 
compact galaxy groups in mocks obtained using three SAMs ap¬ 
plied to the same merger trees derived from the Millennium sim¬ 
ulation. They found that the fraction of compact groups that were 
not dense quartets in real space varied from 24% to 41% depend¬ 
ing on the SAM. In Snaith et al. (2011) four SAMs (two Durham 
and two Munich flavours) all based upon trees extracted from the 
Millennium simulation have been compared with a special focus on 
the luminosity function of galaxy groups - highlighting differences 
amongst models, especially for the magnitude gap distribution be¬ 
tween first- and second-ranked group galaxies. Stripped down ver¬ 
sions of three models (again utilising identical merger trees) have 
been investigated in great detail by De Lucia et al. (2010): they 
studied primarily various assumptions for gas cooling and galaxy 
mergers and found that different assumptions in the modelling of 
galaxy mergers can result in significant differences in the timings of 
mergers, with important consequences for the formation and evolu¬ 
tion of massive galaxies. Most recently Lu et al. (2014) used merger 
trees extracted from the Bolshoi simulation (Klypin et al. 2011) as 
input to three SAMs. They conclude that in spite of the signifi¬ 
cantly different parameterizations for star formation and feedback 
processes, the three models yield qualitatively similar predictions 
for the assembly histories of galaxy stellar mass and star forma¬ 
tion over cosmic time. Note that all three models in this study were 
tuned and calibrated to the same observational data set, i.e. the stel¬ 
lar mass function of local galaxies. Additionally, it should not go 
unmentioned that a lot of effort has gone into comparing SAMs to 
the results of cosmological hydrodynamic simulations, either using 
the same merger trees for both (Yoshida et al. 2002; Helly et al. 
2003; Cattaneo et al. 2007; Saro et al. 2010; Hirschmann et al. 
2012; Monaco et al. 2014) or not (Benson et al. 2001; Lu et al. 
2011) - yielding galaxy populations with similar statistical proper¬ 
ties with discrepancies primarily arising for cooling rates, gas con¬ 
sumption, and star formation efficiencies. For an elaborate review 
of semi-analytical models in relation to hydrodynamic simulations 
please refer to Somerville & Dave (2014). 

In addition to the maturation of SAMs over the past ten years, 
the field of halo models of galaxy clustering has produced other 
powerful techniques for associating dark matter haloes with galax¬ 
ies: the halo occupation distribution (HOD; ling et al. 1998; Cooray 
& Sheth 2002), as well as the complementary models of the con¬ 
ditional luminosity function (Yang et al. 2003) and subhalo abun¬ 
dance matching (Conroy et al. 2006). HOD models have been de¬ 
veloped to reproduce the observed real-space and redshift-space 
clustering statistics of galaxies as well as their luminosity, colour, 
and stellar mass distributions, and they have clarified the simi¬ 
larities and differences between ‘central’ and ‘satellite’ galaxies 
within dark matter haloes. The major difference between SAMs 
and HODs is that the former model physical processes with merger 
trees whereas the latter are relating numerical data (for a given red- 
shift) to observations in a statistical manner, i.e. they bypass an ex¬ 
plicit modelling of the baryonic physics and rely on a statistical de¬ 
scription of the link between dark matter and galaxies. While SAMs 
are therefore guaranteed to return a galaxy population that evolves 
self-consistently across cosmic time, HOD models - by construc¬ 


tion - provide an accurate reproduction of the galactic content of 
haloes. In what follows, we refer to both of them as galaxy forma¬ 
tion models, unless we want to highlight their differences, and in 
those cases we again refer to them as SAM or HOD model. 

In this work - emerging out of the ‘nIFTy cosmology’ work¬ 
shop^ - we are continuing previous comparison efforts, but sub¬ 
stantially extending the set of galaxy formation models: 14 models 
are participating this time, 12 SAMs and 2 HOD models. We are 
further taking a slightly different approach: we fix the underlying 
merger tree and halo catalogue but, for the initial comparison pre¬ 
sented here, we allow each and every model to use its favourite 
parameter set, i.e. we keep the same parameter set as in their ref¬ 
erence model based upon their favourite cosmological simulation 
and merger tree realisation.® By doing this we are deliberately not 
directly testing for different implementations of the same physics, 
but rather we are attempting to gauge the output scatter across mod¬ 
els when given the same cosmological simulation as input. Fur¬ 
ther, this particular simulation and its trees are different from the 
ones for which the models were originally developed and tested. 
By following this strategy, we aim at testing the variations across 
models outside their native environment and without re-tuning. In 
that regard, it also needs to be mentioned that any galaxy forma¬ 
tion model involves a certain level of degeneracy with respect to 
its parameters (Henriques et al. 2009). This is further complicated 
by the fact that different models are likely using different param¬ 
eterizations for the included physical processes. And while it has 
been shown in some of the previously undertaken comparisons that 
there is a certain level of consistency (Fontanot et al. 2009; Lu et al. 
2014), there nevertheless remains scatter. A study including model 
re-calibration will form the next stage of this project and will be 
presented in a future work. 

The remainder of the paper is structured as follows: in Sec¬ 
tion 2 we briefly present the underlying simulation, the halo cata¬ 
logue and the merger tree; we further summarize five different halo 
mass definitions typically used. In Section 3 we give a brief descrip¬ 
tion of galaxy formation models in general; a detailed description 
of each individual model is reserved for Appendix A which also 
serves as a review of the galaxy formation models featured here. 
Some general comments on the layout and strategy of the compari¬ 
son are given in Section 4. The results for the stellar mass function - 
the key property studied here - are put forward in Section 5 whereas 
all other properties are compared in Section 6. We present a discus¬ 
sion of the results in Section 7 and close with our conclusions in 
Section 8. 


2 THE PROVIDED DATA 

The halo catalogues used for this paper are extracted from 62 snap¬ 
shots of a cosmological dark-matter-only simulation undertaken 
using the Gadget-3 W-body code (Springel 2005) with initial 
conditions drawn from the WMAP7 cosmology (Komatsu et al. 
2011, Dm = 0.272, Da = 0.728, Dt = 0.0455, erg = 0.807, 
h = 0.7, Us — 0.96). We use 270® particles in a box of co¬ 
moving width 62.5h~^ Mpe, with a dark-matter particle mass 
of 9.31 X 10®/t“®M©. Haloes were identified with SUBFIND 
(Springel et al. 2001a) applying a Friends-Of-Friends (FOF) link¬ 
ing length of fo = 0.2 for the host haloes. Only (sub-)haloes with at 

^ http :/ /popia . ft . uam.es/nIFTyCosmology 
® Note that not all models entering this comparison were designed ah initio 
to work with N -body trees. 
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least 20 particles were kept. The merger trees were generated with 
the MergerTree code that forms part of the publicly available 
AHF package (Knollmann & Knebe 2009; Gill et al. 2004). For a 
study of the interplay between halo finder and merger tree we refer 
the interested reader to Avila et al. (2014). 

As we run the 14 galaxy formation models with their standard 
calibration, we have to consider five different mass definitions and 
provide them in the halo catalogues; three of them are based upon 
a spherical-overdensity assumption (Press & Schechter 1974) 


Att Q 

Mref(< i?ref) = Aref Pref 77ref , (1) 

and two more on the FOF algorithm (Davis et al. 1985). The 
masses supplied are summarized as follows: 


• M{of. 

• Mhnd' 

• M200c'- 

• M200m- 

• 7Wbn98: 


the mass of all particles inside the FOF group 
the bound mass of the FOF group 
Aref — 200, Pref — Pc 

Aref — 200, Pref — pb 

Aref = AbN 9S, Pref = Pc 


where pc and pb are the critical and background density of the 
Universe, respectively, both of which are functions of redshift and 
cosmology. Abn 98 is the virial factor as given by Eq.(6) in Bryan 
& Norman (1998), and i?ref is the corresponding halo radius for 
which the interior mean density matches the desired value on the 
right-hand side of Eq. (1). 

Note that SUBFIND returns exclusive masses, i.e. the mass of 
host haloes does not include the mass of its subhaloes. Further, the 
aforementioned five mass definitions only apply to host haloes; for 
subhaloes SUBFIND always returns the mass of particles bound 
to it. 


3 THE GALAXY FORMATION MODELS 

The galaxy formation models used in this paper follow one of two 
approaches, which are either ‘semi-analytic model’ (SAM) or ‘halo 
occupation distribution’ (HOD) based models. We give a brief out¬ 
line of both below. A far more detailed description of the semi- 
analytic galaxy formation method can be found in Baugh (2006), 
whereas Cooray & Sheth (2002) provides a review of the formalism 
and applications of the halo-based description of non-linear gravi¬ 
tational clustering forming the basis of HOD models. 

Both model techniques require a dark matter halo catalogue, 
derived from an W-body simulation as described above or produced 
analytically, and take the resultant halo properties as input. Further¬ 
more, SAMs require that the haloes be grouped into merger trees of 
common ancestry across cosmic time. For SAMs, the merger trees 
describing halo evolution directly affect the evolution of galaxies 
that occupy them. HOD models in their standard incarnation do not 
make use of the information of the temporal evolution of haloes; 
HOD models rather provides a mapping between haloes and galax¬ 
ies at a given redshift. However, the utilization of merger trees is 
in principle possible and there are recent advances in that direction 
(Yang et al. 2012; Behroozi et al. 2013). 

The models and their respective reference are summarized in 
Table 1 where we also list their favourite choice for the mass defi¬ 
nition applied for all the results presented in the main body of the 
paper. Some of the models also applied one (or more) of the alter¬ 
native mass definitions and results for the influence on the results 
are presented in Appendix B. The table also includes information 


on whether the models included stellar masses and/or a luminosity 
function during the calibration process. Note that the calibration is 
not necessarily limited to either or both of these galaxy properties. 
And we also chose to provide the assumed initial mass function 
(IMF) in that table which has influences on the stellar masses (and 
star formation rates) presented later in the paper. However, for more 
details about the models please refer to Appendix A. 


3.1 Semi-analytic galaxy formation models 

Semi-analytic models encapsulate the main physical processes gov¬ 
erning galaxy formation and evolution in a set of coupled parame- 
terised differential equations. In these equations, the parameters are 
not arbitrary but set the efficiency of the various physical processes 
being modelled. Given that these processes are often ill-understood 
in detail, the parameters play a dual role of both balancing efficien¬ 
cies and accommodating the (sometimes significant) uncertainties. 
Their exact values should never be taken too literally. However, a 
sensible model is usually one whose parameters are set at an order- 
of-magnitude level of reasonableness for the physics being repre¬ 
sented. All models are calibrated against a key set of observables, 
however exactly which observables are used changes from model 
to model. 

In practice, the semi-analytic coupled equations are what de¬ 
scribe how baryons move between different reservoirs of mass. 
These baryons are treated analytically in the halo merger trees and 
followed through time. The primary reservoirs of baryonic mass 
used in all models in this paper include the hot halo gas, cold 
disk gas, stars, supermassive black holes, and gas ejected from 
the halo. Additionally, different models may more finely delineate 
these mass components of the halo/galaxy system (e.g. breaking 
cold gas into HI and H 2 ) or even add new ones (e.g. the intra-cluster 
stars). 

The physics that a semi-analytic model will try to capture can 
typically be broken into the processes below: 

• Pristine gas infall: As a halo grows its bound mass increases. 
Most SAMs assume that new dark matter also brings with it the 
cosmic baryon fraction of new baryonic mass in the form of pristine 
gas. This gas may undergo heating as it falls onto the halo to form 
a hot halo, or it may sink to the centre along (cold) streams to feed 
the galaxy. At early times, infall may be substantially reduced by 
photoionization heating. 

• Hot halo gas cooling: A hot halo of gas around a galaxy will 
lose energy and the densest gas at the centre will cool and coalesce 
onto the central galaxy in less than a Hubble time. It is usually 
assumed that this cooling gas conserves angular momentum, which 
leads to the formation of a cold galactic disk of gas, within which 
stars can form. 

• Star formation in the cold gas disk: If the surface density of 
cold gas in the disk is high enough, molecular clouds will collapse 
and star forming regions will occur. The observed correlation be¬ 
tween star formation rate density and cold gas density or density 
divided by disk dynamical time can be applied to estimate a star 
formation rate (Kennicutt 1998). 

• Supernova feedback and the production of metals: Very mas¬ 
sive new stars have short lifetimes and will become supernova on 
very short timescales. The assumed initial mass function (IMF) of 
the stars will determine the rate of these events, which will return 
mass and energy back into the disk, or even blow gas entirely out 
into the halo and beyond (Dekel & Silk 1986). Supemovae are also 
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Table 1. Participating galaxy formation models. Alongside the name (first column) and the major reference (last column) we list the applied halo mass 
definition (second column) and highlight whether the calibration of the model included the stellar mass function (SMF) and/or a luminosity function (LF). 
Note that the respective model calibration is not necessarily limited to either of these two quantities. We additionally list the assumed initial mass function 
(IMF). For more details please refer to the individual model descriptions in the Appendix A. 


code name 

mass 

calibration 

IMF 

reference 

Galacticus 

MbN98 

SMF 

LF 

Chabrier 

Benson (2012) 

GalICS-2.0 

Mbnd 

SMF 

— 

Kennicutt 

Cattaneo et al. (in prep.) 

Morgana 

Mfof 

SMF 

— 

Chabrier 

Monaco et al. (2007) 

Sag 

Mbnd 


LF 

Salpeter 

Gargiulo et al. (2014) 

SantaCruz 

MbN98 

SMF 

— 

Chabrier 

Somerville et al. (2008) 

ySAM 

M20OC 

SMF 

— 

Chabrier 

Lee & Yi (2013) 

Durham flavours: 

Galform-gp14 

Mbnd 

— 

LF 

Kennicutt 

Gonzalez-Perez et al. (2014) 

Galform-kb06 

Mbnd 

— 

LF 

Kennicutt 

Bower et al. (2006) 

Galform-kf08 

Mbnd 

— 

LF 

Kennicutt 

Font et al. (2008) 

Munich flavours: 

DLB07 

M2OOC 

— 

LF 

Chabrier 

De Lucia & Blaizot (2007a) 

Lgalaxies 

M20OC 

SMF 

LF 

Chabrier 

Henriques et al. (2013) 

Sage 

M2OOC 

SMF 

— 

Chabrier 

Croton et al. (2006) 

HOD models: 

Mice 

Mfof 


LF 

’diet’ Salpeter 

Carretero et al. (2015) 

Skibba 

M2OOC 

— 

LF 

Chabrier 

Skibba & Sheth (2009) 


the primary channel to produce and move metals in and around the 
galaxy/halo system. 

• Disk instabilities: Massive disks can huckle under their own 
weight. Simple analytic arguments to estimate the stability of a 
disk can be applied to modelled galaxies. Unstable disks will form 
bars that move mass (stars and/or gas) and angular momentum in¬ 
ward. This redistribution of mass facilitates a change in morphol¬ 
ogy where stars pile up in the centre, leading to the growth of a 
bulge or pseudo-bulge (Efstathiou et al. 1982; Mo et al. 1998). 

• Halo mergers: With time haloes will merge, and thus their 
occupant galaxies will as well. Modem cosmological simulations 
readily resolve structures within stmctures, the so-called subhalo 
(and hence satellite) population. However, not all models use this 
information, instead treating the dynamical evolution of substruc¬ 
tures analytically. Subhaloes can either be tidally destroyed or 
eventually merge. It also happens that subhaloes will be lost below 
the mass resolution limit of the simulation, and different models 
treat such occurrences differently. It is common to keep tracking 
the now ‘orphan’ satellites until they merge with the central galaxy. 

• Galaxy mergers and morphological transformation: Galaxy 
mergers are destructive if the mass ratio of the two galaxies is close 
to unity (major merger). In this case the disks are usually destroyed 
to form a spheroid, however, new accreted cooling gas can lead to 
the formation of a new disk of stars around it. Less significant merg¬ 
ers (minor mergers) will simply cause the larger galaxy to consume 
the smaller one. Note that this is a simplified picture and the de¬ 
tailed implementation in a SAM might be based upon more refined 
considerations (e.g. Hopkins et al. 2009) 

• Merger and instability induced starbursts: Mergers and insta¬ 
bilities can also produce significant bursts of star formation on short 
timescales (Mihos & Hernquist 1994b,a, 1996). These are often 
modelled separately from the more quiescent mode of star forma¬ 
tion ongoing in the disk. 

• Supermassive black holes: Most modern SAMs model the for¬ 


mation of supermassive black holes at the centres of galaxies. Such 
black holes typically grow through galaxy merger-induced and/or 
disk instability-induced cold disk gas accretion and/or by the slow 
accretion of hot gas out of the halo. 

• Active galactic nuclei: Accretion of gas onto a supermassive 
black hole will trigger an active phase. For rapid growth this will 
produce a so-called ‘quasar mode’ event, typically occurring from 
galaxy mergers. For more quiescent growth the ‘radio mode’ may 
occur. Feedback resulting from the latter has been used by mod¬ 
ellers to shut down the cooling of gas into massive galaxies, ef¬ 
fectively ceasing their star formation to produce a ‘red and dead’ 
population of ellipticals, as observed. Some models also include 
outflows driven by the ‘quasar mode’. 


It is important to note that not all the processes above are pa- 
rameterised. Often simple analytic theory will provide a reasonable 
approximation for the behaviour of the baryons in that reservoir and 
their movement to a different reservoir. In other cases an observed 
correlation between two properties already predicted by the model 
can be used to predict a new property. However, sometimes there 
is little guidance as to how a physical process should be captured. 
In such cases a power-law or other simple relationship is often ap¬ 
plied. In many cases parameterization of SAM recipes is also based 
upon results from numerical experiments. And as we learn more 
about galaxies and their evolution such prescriptions are updated 
and refined, in the hope that the model overall becomes a better 
representation of the real Universe. 

Fach model used in this paper is described in more detail in 
Appendix A. There the reader can find references that provide a 
complete list of the baryonic reservoirs used and a description of 
the equations employed to move baryons between them. 
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3.2 Halo occupation distribution models 

Halo models of galaxy clustering offer a powerful alternative to 
the semi-analytic method to produce large mock galaxy catalogues. 
Therefore and for comparison to the SAMs, respectively, we in¬ 
clude them in our analysis. As described in the introduction, halo 
models may be of three types: halo occupation distribution (HOD); 
conditional luminosity function (CLF); subhalo abundance match¬ 
ing (SHAM); or some combination or extension of these. Although 
these models are not equivalent, they tend to have similar infer¬ 
ences and predictions for most galaxy properties and distributions 
as a function of halo mass and halocentric distance. The two halo 
models (i.e. MICE and Skibba) in this paper are best described 
as HOD models - even though they both also incorporate subhalo 
properties and hence feature a SHAM component. 

The primary purpose of halo models is to statistically link the 
properties of dark matter haloes (at a given redshift) to galaxy prop¬ 
erties in a relatively simple way. In contrast with SAMs, physi¬ 
cal relations and processes are inferred from halo modeling rather 
than input in the models. By utilizing statistics of observed galaxies 
and simulated dark matter haloes, halo modelers describe the abun¬ 
dances and spatial distributions of cenhal and satellite galaxies in 
haloes as a function of host halo mass, circular velocity, concen¬ 
tration, or other properties, and thus provide a guide for and con¬ 
straints on the formation and evolution of these galaxies. 

Given the halo mass function, halo bias, and halo density pro¬ 
file, HOD models naturally begin with the halo occupation distribu¬ 
tion, P{N\M, c\s, x), where M is the host halo mass, c\s refers to 
central or satellite galaxy status, x refers to some galaxy property 
such as stellar mass or star formation rate. As in SAMs, most mod¬ 
els implicitly or explicitly assume that the cenhal object in a halo is 
special and is (mostly) more massive than its satellites. The mean 
central galaxy HOD, (NcjM), follows an error function, which as¬ 
sumes a lognormal distribution for the central galaxy luminosity (or 
stellar mass) function at fixed halo mass. The mean satellite galaxy 
HOD, (Ns I M ), approximately follows a power-law as a function of 
(M/Mi)“, where the parameter Mi oc Mmin and determines the 
critical mass above which haloes typically host at least one satellite 
within the selection limits. CLF models have different parameteri- 
zations, but the CLFs (or conditional stellar mass functions) may be 
integrated to obtain HODs. All HOD parameters may evolve with 
redshift, though in practice, the stellar mass-halo mass relation, for 
example, evolves very little especially at the massive end. However, 
the occupation number of satellites as a function of mass evolves 
significantly. 

HOD models are constructed to reproduce conditional distri¬ 
butions and clustering of galaxies from their halo mass distribu¬ 
tion, but modeling choices and choices of which constraints to use 
result in different models having different predictions for relations 
between galaxies and haloes. In addition, one must make decisions 
about how and whether to model halo exclusion (the fact that haloes 
have a finite size and cannot overlap so much so that one halo’s 
center lies within the radius of another halo), scale-dependent bias, 
velocity bias, stripping and disruption of satellites, buildup of the 
intracluster light, etc. One can incorporate subhalo properties and 
distributions, as has been done here (which makes the MICE and 
Skibba models HOD/SHAM hybrids), to model orphan satellites 
as well. For bimodal or more complex galaxy properties, such as 
color, star formation rate, and morphology, there is no unique way 
to model their distributions: hence models have different predic¬ 
tions for the quenching and structural evolution of central and satel¬ 
lite galaxies. Models that incorporate merger trees may have differ¬ 


ent predictions for galaxy growth by merging vis-a-vis in situ star 
formation; however, the two models presented here are applied to 
each snapshot individually. 

3.3 Orphan galaxies 

‘Orphan’ galaxies (Springel et al. 2001a; Gao et al. 2004; Guo et al. 
2010; Frenk & White 2012) are galaxies which, for a variety of rea¬ 
sons (such as a mass resolution limit and tidal stripping), no longer 
retain the dark matter halo they formed within. Some of the pro¬ 
cesses that create orphans are physical. For instance, when a galaxy 
and its halo enter a larger host halo tidal forces from the latter will 
lead to stripping of the galaxy’s halo. But there are also numer¬ 
ical issues leading to orphan galaxies: halo finders as applied to 
the simulation data have intrinsic problems identifying subhaloes 
close to the centre of the host halo (see, e.g., Knebe et al. 2011; 
Muldrew et al. 2011; Onions et al. 2012) and therefore the galaxy’s 
halo might disappear from the merger tree (see Srisawat et al. 2013; 
Avila et al. 2014). As, unless they have merged, galaxies do not 
simply disappear from one simulation snapshot to the next, the ma¬ 
jority of the SAMs deal with this by keeping the galaxy alive and 
in their catalogues. They evolve associated properties such as posi¬ 
tion and velocity in different ways: some teams freeze these at the 
values they had when the host halo was last present while others 
integrate their orbits analytically or tag a background dark matter 
particle and follow that instead. Note that information about the 
motion of dark matter particles was not supplied to the modellers 
during the comparison and so this latter method was not available 
in this study. For that reason we do not include any analyses which 
are affected by the choice of how to assign positions and velocities 
to orphans. 

3.4 Calibration 

For the initial comparison presented here we do not require a com¬ 
mon calibration, but rather each semi-analytical model is used ‘as 
is’. In particular this means that each model has been run with its 
preferred physical prescriptions and corresponding set of parame¬ 
ters, as detailed in the Appendix, without specific retuning. One of 
the purposes in doing this is to show the importance of providing 
adequate calibraton when applying the models. We show in Sec¬ 
tion 5 that the raw scatter in the stellar mass function is very large 
but reduces considerably when models use their preferred merger 
tree/halo mass definition and initial-mass-function. A detailed com¬ 
parison of models when retuned to match the nIFTy merger trees 
and halo catalogues will be the topic of a future paper. 


4 THE COMPARISON 

First we present our methodology and establish the terminology 
used throughout the remainder of the paper. We illustrate in Fig. 1 
how galaxies and the haloes they live within are connected. A cir¬ 
cle represents a dark matter halo, a horizontal line a galaxy. A solid 
arrow points to the host halo of the galaxy whereas a dashed ar¬ 
row points to the main halo the galaxy orbits within. Note that any 
substructure hierarchy has been flattened to one single level, i.e. 
sub-subhaloes will point to the highest level main halo. The sketch 
also includes orphan galaxies which will have a pointer to its ‘last’ 
dark matter host halo. 

To facilitate comparison, in all subsequent plots that make use 
of halo mass we (arbitrarily) chose Mbnd regardless of the mass 
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Figure 1. Sketch illustrating the population of galaxies residing inside a 
dark matter ‘main halo’ (outermost circle). Each galaxy can have its own 
dark matter halo which we will refer to as ‘host halo’ (or simply ‘halo’). 
Galaxies can also be devoid of such a ‘galaxy host halo’ and then they are 
tagged as an ‘orphan’ galaxy. 


definition used internally by each model. Each galaxy in the sup¬ 
plied galaxy catalogues points to its host halo (either the present 
one, if it still exists, or its last one in the case of orphans) and 
this link is used to select Mbnd from the original halo catalogue 
- irrespective of the mass definition used in the actual model. This 
allows us to directly compare properties over a single set of halo 
masses without concern that the underlying dark matter framework 
has shifted slightly or is evolving differently with redshift. 

We have further limited all of the comparisons presented here 
to galaxies with stellar masses M* > - a mass thresh¬ 

old appropriate for simulations with a resolution comparable to the 
Millennium simulation (see Guo et al. 2011). Additionally, when 
connecting galaxies to their haloes we have applied a mass thresh¬ 
old for the latter of Mhaio > corresponding to ap¬ 

proximately 100 dark matter particles for the simulation used here. 

When interpreting the plots we need to bear in mind - as can 
be verified in Appendix A - that a great variety of model cali¬ 
brations exist: some models tune to stellar mass functions (SMF), 
some models to luminosity functions (LF), and the HOD models 
to the clustering properties of galaxies. To facilitate the differentia¬ 
tion between calibrations in all subsequent plots we chose colours 
for the models as follows: models calibrated using only SMFs are 
presented in blue, models calibrated using LFs in red, and models 
using a combination of both in green; the 2 HOD models are shown 
in black, but with different linestyles. 

Note that not all models necessarily appear in all plots. For in¬ 
stance, the Skibba hod model only provided a galaxy catalogue 
for redshift 2 = 0 and hence is not shown in evolutionary plots. 
Some models do not feature orphans and so do not appear in the 
corresponding plots. 

For the comparison presented in this first paper we now focus 
on the stellar mass function (SMF) to be presented in the following 
Section 5. There we will also investigate several origins for model- 
to-model variations seen not only for the SMF, but also for other 
galaxy properties to be presented in Section 6, e.g. star formation 
rates, number density of galaxies, orphan fractions, and the halo 
occupation. We deliberately exclude luminosity-based properties as 




Figure 2. Stellar mass function at redshift 2 = 0 (top) and 2 = 2 (bot¬ 
tom). Each model used its prefeired mass definition and initial stellar mass 
function. 

they introduce another layer of modelling, i.e. the employed stellar 
population synthesis and dust model. 


5 STELLAR MASS LUNCTION 

A key property, used by many of the models presented here to con¬ 
strain their parameters, is the stellar mass function of galaxies: we 
therefore dedicate a full section to its presentation and discussion. 
It is shown in Fig. 2 both at redshift 2 = 0 (top panel) and 2 = 2 
(bottom panel). This plot indicates that there is quite a range in both 
galaxy abundance and mass (influenced by star formation rate and 
star formation history) across the models. These differences are ap¬ 
parent at 2 = 0 where the models are calibrated, but are even more 
pronounced at redshift 2 = 2. At 2 = 0 the model results vary 
in amplitude by around a factor of 3 in the main and exhibit high 
mass cut-offs of varying steepness and position; at 2 = 2 the dif¬ 
ferences in amplitude are even larger, reflecting a broad variation 
in the location of the peak in star formation rates (presented in the 
next section). 

The HOD model MICE lies within the range of stellar mass 
functions provided by the SAMs as does the SKIBBA model above 
10®'®/i“^M©. At 2 = 2 the SKIBBA model does not provide a 
return while the MICE model features amongst the models with 
the largest number of high mass galaxies, i.e. the MORGANA and 
Sag models. While in MORGANA the overproduction of massive 
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Figure 3. Stellar mass function at redshift z = 0 for models that (also) 
returned galaxy catalogues using M 200 C as the mass definition. To be com¬ 
pared against the upper panel of Fig. 2 


Figure 4. Stellar mass function at redshift z = 0 after applying a correction 
for the applied IMF, i.e. models have been coirected towards a Chabrier 
IMF. To be compared against the upper panel of Fig. 2 


galaxies is connected to the inefficiency of the chosen AGN feed¬ 
back implementation to quench cooling in massive haloes, in SAG 
and Mice it could additionally be related to the assumption of 
a Salpeter IMF which implies a higher mass estimate than for 
a Chabrier IMF (see Section 5.2 below). The mass function for 
Lgalaxies at z = 2 is lower than any other model due to the 
delayed reincorporation of gas ejected from supernova feedback 
that shifts star formation in low mass galaxies to later times (Hen- 
riques et al. 2013). At both redshifts the Galacticus model dis¬ 
plays a bump in the stellar mass function around due 

to the matching of feedback from active galactic nuclei and super¬ 
novae. For completeness we also checked that the scatter seen here 
basically remains unchanged when restricting the analysis to (non- 
)central galaxies and (non-)orphans, respectively. 

The differences seen here are huge, especially at the high- 
mass end, even when models have implemented the same physical 
phenomena such as supernova and AGN feedback. For instance, 
Lgalaxies and Galacticus both allow the black hole to ac¬ 
crete from the hot halo, with associated jets and bubbles producing 
‘radio mode’ feedback: however, the mass of the largest galaxies 
differs by around an order of magnitude at redshift 2 ; = 0. In order 
to understand how much of this difference arises from the different 
physical implementations, we first need to consider other factors 
that may influence the results. For example, the models: 

a) use a variety of halo mass definitions; 

b) use different initial-mass functions (IMFs); 

c) have been taken out of their native environment, i.e. they have 
been applied to a halo catalogue and tree structure that they were 
not developed or tested for; 

d) have not been re-calibrated to this new setup; and 

e) have not been tuned to the same observational data. 

In the following sub-sections we will address points a-c) in more 
detail. Points d) and e) are more complex and will be left for a 
future study. 

5.1 Mass Definition 

It can be seen from Table 1 that the models participating in this 
comparison applied a variety of different mass definitions (which 
were introduced in Section 2) to define the dark matter haloes that 



Figure 5. Stellar mass function at redshift 2 = 0 for all the models when 
created in the native environment used during the model calibration. All 
curves have been corrected for the IMF towards Chabrier. This figure uses 
the same scale as (and should be compared to) Fig. 4. 

formed their halo merger tree. But as several of the code represen¬ 
tatives also returned galaxy catalogues using mass definitions other 
than their default one, we are able to prepare a plot that shows the 
stellar mass function for M 200 C, i-e- the mass definition for which 
the maximum number of models exist. We show that plot as Fig. 3 
where we see that the effect of changing the mass definition is 
smaller than the model-to-model variation and hence not the pri¬ 
mary source of it. 

Appendix B provides a direct comparison of models for two 
different mass definitions (their standard one and M 2 ooc)- That ap¬ 
pendix further shows its influence on other galaxy properties such 
as the stellar-to-halo mass ratio and the number and star formation 
density evolutions. 

5.2 IMF Correction 

An additional source of scatter is that the models assumed vari¬ 
ous initial stellar mass functions. Hence we transformed the stellar 
masses returned by each model to a unified Chabrier IMF (Chabrier 
2003). For that we used the following equations (Mitchell et al. 
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Figure 6. Comparing the stellar mass functions for each model as given in its native environment (thick lines) and when applied with the same parameters to 
the nIFTy data (thin lines). 


2013; Bell &de Jong 2001): 


-0.240 

= logio(A/f"*-®"'P"*®") -0.090 

logio(Mf = logio(Mf™"““) +0.089 

( 2 ) 

Note that this is only a rough correction, as these numbers depend 
on the stellar population synthesis (SPS) model, age and metallicity 
of the simple stellar population and on looking to one or several 
bands when estimating stellar masses from broad band photometry. 
The models have been corrected as follows: 

• GalICS- 2.0: tuned to observations w/ Chabrier IMF; 

• Galform: Kennicut —>■ Chabrier; 

• Mice: diet-Salpeter —?■ Chabrier; 

• Sag: Salpeter —>■ Chabrier; 

noting that we left GalICS- 2.0 untouched because this model 
tuned its parameters to an observational data that itself assumed 
already a Chabrier IMF. 

In Fig. 4 we show the resulting stellar mass function for all 
models where we notice again that the scatter is only slightly 
reduced. Some additional information is again provided in Ap¬ 
pendix C. 


5.3 Model Environment 

While the whole idea of the comparison presented in this paper is 
to apply galaxy formation models to the same halo catalogues and 
merger trees coming from a unique cosmological simulation, we 
have seen that there is a non-negligible scatter across properties. 
This scatter is larger than in previous comparison projects which 
encompassed fewer models (e.g. Kimm et al. 2009; Fontanot et al. 
2009, 2012; Contreras et al. 2013; Lu et al. 2014). As we could 
neither attribute the increased variations to halo mass definitions 
nor the assumed initial stellar mass functions we are now going to 
show that this in part comes from taking models out of their native 
environment. The majority of the models have been designed using 
a certain simulation and tree structure. However for the compar¬ 
ison presented here this environment has often been substantially 
changed, and model parameters have not been adjusted to reflect 
this (as mentioned before, such a recalibration will form part of the 
follow-up project). The effect of this approach can be appreciated 
in the two following plots. 

Fig. 5 shows the stellar mass functions at redshift 2 = 0 for all 
models as given when applied to the simulation and merger trees 
used for calibration', we refer to this setup as ‘native environment’. 
These data points were directly provided by the code representa¬ 
tives, then converted to a common IMF, as in the previous section, 
and plotted on the same axes/scale as for Fig. 2. The agreement be¬ 
tween the models is now much improved indicating that the main 
part of the scatter seen before is due to models being applied to 
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simulation data they were not adapted to (and that might even fea¬ 
ture a different cosmology). There are still a few outliers on the 
low-mass end, e.g. GalICS- 2.0: this model was calibrated on the 
Ultra VISTA SMFs at different redshifts, not on local mass func¬ 
tions; and models calibrated on high -2 data tend to underestimate 
the low-mass end of the local SMF. 

To underline the influence of the merger trees to the model 
results we directly compare in Fig. 6 the native (thick lines) to the 
nIFTy (thin lines) stellar mass function where each panel represents 
one model. This plot quantifies the sensitivity of the model to the 
underlying simulation and merger tree. It shows that recalibration 
is required whenever a new simulation is to be used. While this 
is common practice within the community, its necessity has been 
shown here for the first time. A forthcoming companion paper will 
further address the influence of the applied observational data set 
to the remaining model-to-model scatter. 


6 GALAXIES AND THEIR HALOES 

In this section we extend the comparison to several additional prop¬ 
erties including star formation rate, the stellar mass fraction, the 
number density (evolution) of galaxies, and the relation between 
galaxies and their dark matter haloes. 

6.1 Star formation rate 

The stellar mass of a galaxy studied in the previous section depends 
upon the evolution of its star formation rate (SFR). Therefore we 
now turn to the history of the star formation rate across all consid¬ 
ered models. In Fig. 7 we show the redshift evolution, noting that 
all the curves in this plot have been normalized by their redshift 
2 = 0 values (which are given in the third column of Table 2). 
In this way we separate trends from absolute differences. An un¬ 
normalized version of Fig. 7 can be found in Appendix D. Remem¬ 
ber that the HOD model Skibba does not produce high redshift 
outputs and so appears neither in Fig. 7 nor in the SFR columns of 
the accompanying Table 2. 

For the SAMs the peak of star formation is about redshift 
2 ~ 2 — 3 followed by a rapid decrease at late times - in agree¬ 
ment with observations and the uncertainties seen within them (e.g. 
Madau & Dickinson 2014). But amongst these models there are 
also differences: in Lgalaxies, for instance, the peak is at smaller 
redshifts while for Galacticus it is at earlier times; and the HOD 
model Mice shows a relatively high SFR at low redshifts, i.e. 
Mice stars are formed preferentially later than in the other models. 
These differences are reflected in the fact that, given the redshift 
2 = 0 normalisation in the plot, there are differences in amplitude 
of an order of magnitude at redshift 2 > 6 . 

While the previous figure has shown the integrated star forma¬ 
tion rate, we inspect its redshift 2 = 0 properties more closely in 
Fig. 8 where we present the star-formation-rate distribution func¬ 
tion, i.e. the number density of galaxies in a certain SFR interval. 
From that we see that all models have a similar functional form, but 
that the normalisation of the SFR shows differences of up to a factor 
of three between models. This is reflected in Table 2 where we list 
the total stellar mass formed (second column), the present-day star 
formation rate (third column) and specific star formation rate (last 
column, i.e. the ratio between SFR and total M*). We see that, for 
instance, the GALACTICUS model produced more than three times 
as many stars as Lgalaxies. 

The two previous figures showed the overall star formation 


Table 2. Total stellar mass, star formation rate, and global specific star for¬ 
mation rate in galaxies with M* > lO®h“^M 0 at redshift 2 = 0 (com¬ 
puted as total stellar mass divided by total SFR). 


code name 

0 

O 

SFR 

[lO‘‘h-iM 0 yr-l] 

sSFR 

Gyr-i 

Galacticus 

2.91 

2.22 

0.0761 

GalICS-2.0 

2.73 

0.88 

0.0321 

Morgana 

1.96 

1.21 

0.0614 

Sag 

2.37 

1.15 

0.0486 

SantaCruz 

1.11 

0.53 

0.0475 

ySAM 

1.14 

0.85 

0.0749 

Durham flavours: 

Galform-gp14 

0.98 

0.50 

0.0511 

Galform-kb06 

1.16 

0.51 

0.0442 

Galform-kf08 

1.06 

0.52 

0.0491 

Munich flavours: 

DLB07 

1.76 

0.99 

0.0563 

Lgalaxies 

0.87 

1.07 

0.1234 

Sage 

1.01 

0.83 

0.0815 

HOD models: 

Mice 

1.77 

0.96 

0.0543 

Skibba 

1.49 

n/a 

n/a 



Figure 7. Star formation rate density for galaxies with M* > 10®/r 

as a function of redshift (normalized to the redshift 2 = 0 values listed in 

Table 2). 

rate, but now we focus in Fig. 9 on the specific star formation rate 
(sSFR) as a function of stellar mass M* at redshift 2 = 0. The 
left panel shows the sSFR excluding passive galaxies, i.e. galaxies 
that are not considered ‘star-forming’, which we define as those 
with sSFR < 0.01 Gyr“^. Points shown are the mean values in 
the bin, both for the y- and a;-axes, (which explains why they do 
not start exactly at our mass threshold of IO^/i^^Mq).^ Instead 
of error bars, the right panel of Fig. 9 shows the distribution of 
sSFR values for a mass bin M* € [10^°, lO^^]/i“^M 0 with our 

* Although not shown we also reproduced the plot using medians and 25 
and 75 percentiles, but as these give very similar results we decided to adopt 
mean values for this and all subsequent plots. 
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Figure 9. Specific star formation rate of star-forming galaxies at redshift z = 0. The left panel shows the star formation rate per stellar mass as a function of 
stellar mass M*; points represent mean values binned in both the y and x direction for the star forming sequence of galaxies. The right panel seiwes as a proxy 
for the (omitted) en'or bars: it shows the distribution of specific star formation rates for galaxies in the mass range M* S [10^*^, (indicated by 

the two vertical lines); the vertical dashed line indicates our choice for the passive galaxy fraction threshold. 



Figure 8. Star formation rate distribution function at redshift 2 = 0. 


choice for the passive threshold shown as a vertical dotted line. 
We are aware that the choice of this mass bin for the right panel 
encompasses the ‘knee’ of the stellar mass function, but this right 
panel nevertheless shows that our passive threshold value cuts the 
wing to the left at approximately the same height as the right side. 
Please note that the right panel does not substantially change when 
considering a different mass range, although this is not explicitly 
shown here. 

Fig. 9 reflects what has already been seen in Fig. 7, i.e. there 
is a great diversity in star formation rates across the models irre¬ 
spective of the stellar mass of the galaxy. Bearing in mind the dif¬ 
ferences in stellar mass functions visible here again on the x-axis, 
the curves vary by a factor of about 3 at essentially all masses, with 
the primary difference being in overall normalisation. 

We would like to remark on the interplay between Fig. 7, Ta¬ 
ble 2 and Fig. 9 as at first sight the results seem to be counterintu¬ 
itive. For instance, Galacticus has a much higher star formation 
rate (at all times) than Lgalaxies, yet the specific star formation 
rate is higher for Lgalaxies. This is readily explained by the pres¬ 
ence of, on average, more massive galaxies in GALACTICUS, which 
is confirmed by the stellar mass function presented in Fig. 2. One 
should also bear in mind that the specific star formation rate sSFR 


could be considered a proxy for the (inverse of the) age of a galaxy. 
Therefore, Fig. 9 indicates, for example, that it took galaxies in 
Lgalaxies less time to assemble their stellar mass than galaxies 
in Galacticus. 

We close our discussion of Fig. 9 with the remark that it does 
not change when considering only centrals: the differences across 
models remain unaffected by restricting the analysis to this galaxy 
population. However, the specific star formation rate rises, with ap¬ 
proximately constant ratios between the curves, when moving to 
higher redshifts. 

6.2 Stellar mass fractions 

In Fig. 10 we show the stellar-to-halo mass ratio M,/Mhaio as a 
function of galaxy host halo mass Mhaio for non-orphan galaxies.® 
The layout of this figure is similar to the previous one, i.e. the left 
panel shows the actual mean of the ratios (omitting error bars), 
with the distribution of the values in a galaxy host halo mass bin 
Mhaio e [ 10 ^^, the right serving as a proxy for the 

missing error bars. All models have a similar maximum stellar mass 
fraction of about 0 . 1 , but there is a large spread in the modal value 
and this leads to a difference in overall normalisation of approxi¬ 
mately a factor of 3 - where some of this variation can be attributed 
to the different mass definitions applied (see Appendix B). Note 
that the curves remain unaffected by restricting the data to central 
galaxies only. Further, the distributions shown in the right panel are 
not influenced by the halo mass bin. 

6.3 Number density 

Turning to the galaxies themselves we show in Fig. 11 the evolu¬ 
tion of the number density of galaxies (with stellar mass in excess 
of M, > lO®/i“^M 0 ) as a function of redshift - normalized to 
the number of galaxies at redshift 2 = 0 (provided in Table 3).® It 
is noteworthy that the various models display different evolution¬ 
ary trends of galaxy density. For instance, DLB07 starts with the 

® We omit orphan galaxies as they lack an associated halo mass. 

® Note that as the Skibba model solely provides 2 = 0 data it does not 
appear in the figure. 
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Figure 10. Stellar-to-halo mass ratio for all (non-orphan) galaxies at redshift z = 0. The left panel shows mean values (again in both directions) as a function 
of galaxy host halo mass Mnaio whereas the right panel indicates the (omitted) error bars: it shows the distribution of AT*/Mijaio for galaxy halo masses in 

the range Mjjaio S [ 10 ^^, lO^^]h“^M 0 . 



Figure 11. The number of all galaxies with stellar mass M* > 
10®h“^MQ (normalized to redshift z = 0 values as listed in Table 3) as 
a function of redshift. The inset panel shows a zoom (using a linear j/-axis) 
into the range z S [ 0 , 1 ]. 

largest fraction of galaxies at high redshift whereas GalICS- 2.0 
begins with the lowest fraction - with the difference being more 
than one order of magnitude at redshift z = 6 between these two 
models. We further note (in the inset panel) that some of the models 
- in particular SAGE - have a flat or falling galaxy number density 
between a redshift z = 1 and the present day - this is perfectly al¬ 
lowable as galaxy merging can reduce the galaxy density. Although 
not explicitly shown here, the main features of the plot remain the 
same when restricting the analysis to central galaxies only. 

Fig. 11 should be viewed together with Table 3 as the former 
provides the trend whereas the latter quantifies the normalization 
(at redshift 2 = 0 ); for a combination of both, i.e. an un-normalized 
version of Fig. 11, we refer the reader to Fig. D2 in the Appendix. 
The total number of galaxies with M* > ranges from 

~ 7500 for the Lgalaxies model to «15000 in the DLB07 
model. All models (apart from Skibba) populate all dark matter 
(sub-(haloes found in the simulation down to at least Mhaio ~ 
lO^^/i“^M 0 . Flence, any differences seen here originate from 
lower mass objects. This is confirmed by re-calculating Wcentrai 
applying a halo mass threshold of Mhaio > 2 x lO^^fi“^M 0 (in- 


Table 3. Number of galaxies at redshift 2 = 0 with a stellar mass in excess 
of M* > lO®/r“'^M 0 . (For the MORGANA and SantaCruz models the 
number of oiphans is in fact the number of satellite galaxies, see text.) 


code name 

Agal 

-^central 

-^non —orphan 

-^orphan 

Galacticus 

14255 

7825 

10019 

4236 

GalICS-2.0 

9310 

7462 

9310 

0 

Morgana 

10008 

6186 

6186 

3822 

Sag 

19516 

13571 

16256 

3260 

SantaCruz 

8901 

6682 

6682 

2219 

ySAM 

11138 

7423 

9458 

1680 

Durham flavours: 

Galform-gp14 

8824 

5097 

6098 

2726 

Galform-kb06 

11563 

6669 

7897 

3666 

Galform-kf08 

12116 

6430 

7664 

4452 

Munich flavours: 

DLB07 

15132 

9420 

11897 

3235 

Lgalaxies 

7499 

4792 

6287 

1212 

Sage 

8437 

6588 

8437 

0 

HOD models: 

Mice 

12191 

7286 

10106 

2085 

Skibba 

9203 

5088 

7973 

1230 


stead of the galaxy stellar mass threshold of M* > lO®/i“^M 0 ). 
This process results in 3774 galaxies for all models - a number 
identical to the number of host haloes in the SUBFIND catalogue 
above this mass limit. 

In Table 3 we further divide the galaxies into different pop¬ 
ulations, i.e. centrals, non-orphans, and orphans. The fraction of 
orphan galaxies also shows a spread from a mere 13 per cent for 
Skibba to nearly 37 per cent for the Galform-KFOS model. 
Note that GalICS- 2.0, and SAGE do not feature orphans at all, 
whereas the MORGANA and SantaCruz models, as previously 
mentioned, do not make use of the A-body information for sub¬ 
haloes and hence tag satellite galaxies as orphans - therefore all 
satellite galaxies in these models are technically orphans as only 
central galaxies retain information about their host halo; naturally. 
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-^central = A^non-orphan for these two models. To further explore 
the differences in the abundance of orphan galaxies between mod¬ 
els we show in the upper (lower) panel of Fig. 12 the number (stel¬ 
lar mass) fraction of all galaxies that are classified as non-orphans 
with stellar mass in excess of M* > orbiting inside a 

main halo of given mass Mhaio-^ We can observe some bimodal¬ 
ity here: the two HOD models have the lowest fraction of orphans 
in high-mass main haloes (cf. Table 3) whereas for all other mod¬ 
els orphans form the dominant population, making up between 40 
and 75 per cent of all galaxies at z = 0 within main haloes above 
Afhaio > This trend is also true at higher redshifts 

although we do not explicitly show it here. For those models which 
feature orphans, the variation in the number of orphans is due to 
the various methods of dealing with their eventual fate: over some 
timescale orphans are expected to suffer from dynamical friction 
and merge into the central galaxy of the halo. This timescale can be 
very long in some models (see e.g. De Lucia et al. 2010). The ba¬ 
sic features of the plot do not change when applying a more strict 
threshold for M* as can be verified in the lower panel of Fig. 12 
where we show the stellar mass weighted fraction of non-orphan 
satellite galaxies. However, now the model variations are reduced. 
The difference from unity is the fraction of stellar mass locked up 
in orphans, still as large as 60 per cent for the Galform mod¬ 
els, yet substantially smaller than the number fraction presented in 
the top panel. These differences might be ascribed to the different 
treatment of merger times again. We also observe that for the HOD 
models the orphan contribution has nearly vanished when weighing 
it by stellar mass. 

6.4 Galaxy-Halo Connection 

Lastly, we now turn to how the supplied tree hierarchy has been 
populated with galaxies by each model. To this extent the upper 
panel of Fig. 13 compares the supplied halo mass function (crosses, 
all identified objects down to 100 particles including subhaloes) 
overplotted by the halo mass function as derived from the galaxy 
catalogues returned by the models. Note that the MORGANA and 
SantaCruz models have been omitted due to their treatment of 
subhaloes. While the logarithmic scale of the upper panel masks 
any differences, the lower panel - showing the fractional difference 
of each model’s halo mass function with respect to the supplied 
input mass function ~ indicates that nearly every (sub-)halo found 
in the simulation contains a galaxy. The only exception to this is 
the Skibba model, which has a high incompleteness threshold that 
frequently leaves small haloes empty. For the other models differ¬ 
ences are all below 2 per cent. 

In Fig. 14 we relate the number of galaxies with stellar mass 
M* > to the mass of the dark matter main halo 

they orbit within. We normalize the average number of galaxies 
by the mass of the main halo, i.e. presenting the ‘specific fre¬ 
quency of galaxies’. The solid line running from the upper left 
to the lower right across the plot indicates a frequency of ‘one 
galaxy per main halo’. Haloes to the right of this line essentially 
always contain at least one galaxy while the values to the left are 
indicating the fraction of haloes devoid of galaxies and are chiefly 
driven by incompleteness - and in retrospect justifying our thresh¬ 
old of Mhaio > lO'^^/i“^M 0 for the previous plots. At high halo 
mass, Mhaio > lO^^/i~^M 0 , the specific frequency of galaxies 

As neither GalICS- 2.0, MORGANA, SAGE, nor SantaCruz feature 
orphans, they have been omitted from the plot. 




Figure 12. Number fraction (top panel) and stellar mass fraction (bottom 
panel) of non-orphan galaxies with stellar' mass M* > lO®fi“^M 0 as a 
function of main halo mass Mhaio at redshift z = 0. Only models that 
actually feature orphans are shown here. 


is roughly constant for all models although the occupation num¬ 
ber of haloes varies between them by around an order of magni¬ 
tude. While the upper panel of Fig. 14 shows all galaxies - includ¬ 
ing orphans - the lower panel only shows non-orphan galaxies; we 
clearly see that the differences between models are primarily due 
to the (treatment of) orphan galaxies, though significant differences 
remain even for non-orphans. The difference between the two pan¬ 
els also indicates that orphans are favourably found in higher mass 
haloes while lower mass objects are practically devoid of them, as 
already seen in Fig. 12. 

The halo occupation distributions (as well as clustering prop¬ 
erties) of all the models will be further analyzed in a spin-off project 
of this collaboration (Pujol et al., in prep.). 


7 SUMMARY & DISCUSSION 

We have brought together 14 models for galaxy formation in sim¬ 
ulations of cosmic structure formation, i.e. 12 SAM and 2 HOD 
models. In this inaugural paper we presented the models and under¬ 
took the first comparison where the models applied their published 
parameters (without any recalibration) to the same small cosmolog¬ 
ical volume of (62.5/i“^M0)^ with halo merger trees constructed 
using a single halo finder and tree building algorithm. Hence the 
framework that underpins this study was designed to be the same 
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Figure 13. Mass function of all haloes at redshift 2 = 0 as given in the 
input halo catalogue (crosses, all identified objects down to 20 particles in¬ 
cluding subhaloes) and as recovered from the non-orphan galaxy catalogues 
of each model. The upper panel shows the supplied mass function whereas 
the lower panel shows the fractional difference with respect to the input 
halo catalogue. The upper panel also gives the translation of Muaio to the 
number of particles in the halo as additional a;-axis at the top. 

for each model. This approach allowed us to directly compare the 
galaxy formation models themselves leaving aside concerns about 
cosmic variance, the influence of the halo finder (Avila et al. 2014) 
or tree construction method (Srisawat et al. 2013; Lee et al. 2014). 
However, some teams had to slightly alter this framework in order 
to make the catalogues compatible with the assumptions in their 
methods (see Appendix A). 

All of the contributing teams have been provided with a stan¬ 
dardized dark matter halo catalogue and merger tree; they were 
asked to undertake their currently favoured model and were ex¬ 
plicitly told the underlying cosmology and mass resolution of the 
simulation to be used. They supplied returns in a specified format 
and the analysis was performed on these files using a single com¬ 
mon analysis pipeline. This approach has proved highly successful 
for other related scientific issues in the past (e.g. Frenk et al. 1999; 
Heitmann et al. 2008; Knebe et al. 2011; Onions et al. 2012; Sri¬ 
sawat et al. 2013; Kim et al. 2014). 

This paper should be viewed as the first in a series emerg¬ 
ing out of the nIFTy cosmology workshop®. A number of spin-off 
projects were also initiated at the meeting, including more detailed 

® http://www.popia.ft.uam.es/nIFTyCosmology 




Figure 14. Number of galaxies per halo mass (i.e. ‘specific frequency 
of galaxies’) for galaxies more massive than M* > as a func¬ 

tion of halo mass ATjjaio at redshift 2 = 0. Points plotted are the mean 
values in the bin with respect to both axes and error bars are Itr. The solid 
line running from upper left to lower right represents one galaxy per halo. 
The upper panel shows all galaxies whereas the lower panel focuses on non- 
oiphan galaxies. 


studies of cold vs. hot gas properties, correlation functions, dust 
effects, disk instabilities, and - last but not least - how to define 
a common calibration framework. The results will be presented in 
future papers. 

For the present paper each team applied their published cali¬ 
bration values to the supplied cosmological model as specified in 
Appendix A which also describes each model’s specific choice of 
parameters. Each team tends to use their own personal preference 
of which observables to tune their model to, and these observables 
often require additional processing to produce from the more phys¬ 
ically fundamental quantities studied here. For instance, the calcu¬ 
lation of luminosities requires the adoption of a particular stellar 
population synthesis model as well as a certain dust model. Even 
the derivation of stellar mass demands a choice for the stellar initial 
mass function. We deliberately deferred from studying magnitudes 
to avoid the accompanying layer of complexity; and we left the 
preference for the initial stellar mass function at the modeller’s dis¬ 
cretion. This choice of calibration freedom was made deliberately 
in order not to favour any particular model if we happened to make 
similar (somewhat arbitrary) post-processing choices. 

We have explicitly chosen not to overlay our figures with ob¬ 
servational data for precisely the same reason: such data requires 
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somewhat arbitrary (reverse) conversion from the observed quanti¬ 
ties and this conversion may bias the reader in favour of a partic¬ 
ular model that happened to convert using the same approach (or 
happened to tune their model to this particular observational quan¬ 
tity). We reserve such comparisons to future work where we will 
consider the full range of such conversions and include a careful 
review of the observational literature on this point. But we never¬ 
theless like to remind the reader that each model has been compared 
with a variety of observational data and that these comparisons are 
published. 

Given the variety in models and calibrations, the agreement 
found here is gratifying although a number of discrepancies exist, 
as summarized here (and discussed below). 

Stellar Component The stellar mass functions at 2 = 0 of all the 
models lie within a range of around a factor of 3 in amplitude across 
the faint end of the curve and have somewhat different effective 
breaks at the high mass end. The star formation rate density and 
specific star formation rate by halo mass is broadly similar for all 
the SAMs, with the main difference being in normalisation, which 
can vary by a little under an order of magnitude near the peak of the 
star formation rate density curve and a factor of three elsewhere. 

Galaxies For most of the models considered here, galaxies without 
a surviving dark matter halo (so-called ‘orphan galaxies’) domi¬ 
nate the number counts within each host halo, accounting for be¬ 
tween 40 and 75 per cent of all galaxies in main haloes above 
at redshift 2 = 0. The treatment of these galaxies and 
their eventual fate differs dramatically between the various models 
and is also expected to be strongly dependent upon the resolution 
of the simulation. 

Galaxy-Halo Connection All the models populate the supplied 
trees adequately, i.e. all haloes found in the simulation contain a 
galaxy. We further found that the specific frequency of galaxies, 
i.e. the number of galaxies per halo mass, is constant above the 
completeness limit of the simulation although the average number 
of galaxies per dark matter halo mass varies by around an order 
of magnitude across the models, if orphans are included; otherwise 
the variation is reduced to a factor of about two. 

When interpreting the results one needs to always bear in mind 
that all models were used as originally tuned in the respective ref¬ 
erence paper, i.e. the way this comparison has been designed might 
lead to scatter across models that is larger than the scatter due to 
different implementations of the same physics within them. The 
factors entering into the model-to-model variations seen here are 
differences due to a) models not being tuned to the same observa¬ 
tional constraints, b) models being tuned to different cosmologies, 
c) the choice for the halo mass definition, d) the choice for the ap¬ 
plied initial stellar mass function, and e) models not being opti¬ 
mally tuned (for the merger tree structure at hand). Elaborating on 
these points: 

a) observational constraints: Using different observational con¬ 
straints should not be a primary source of variation, at least as 
long as constraints for all the relevant physical quantities are in¬ 
cluded. For instance, the work of Henriques et al. (2013, using 
SMF-l-Bband-l-Kband constraints from 2 = 3 to 2 = 0) and Hen¬ 
riques et al. (2014, using SMF-^red fraction constraints) lead to 


convergent results, if a proper assessment of the observational un¬ 
certainties is performed; the authors also state that they need the 
combination of properties to arrive at converged likelihood regions 
in parameter space. However, the models included here show an 
even greater variety of (potentially mutually exclusive) observa¬ 
tional constraints and hence we cannot exclude that those differ¬ 
ences contribute significantly to the scatter. 

b) cosmology: Differences due to cosmology can often (but not 
always) be absorbed by re-tuning the physical parameters of the 
model. However, without re-tuning (as here), cosmology can make 
a big difference; this can be seen, for instance, from the significant 
changes in parameter values required to get the same stellar mass 
function at 2 = 0 for different cosmologies (see Wang et al. 2008; 
Guo et al. 2013; Gonzalez-Perez et al. 2014). 

c) halo mass definition: As can be verified in Table 1 another dif¬ 
ference across models is the applied definition for halo masses. Fur¬ 
thermore these differences can increase with redshift: for instance, 
M 20 QC only depends on the evolution of pc{z) whereas Mbn 9 s has 
an additional dependency on cosmology encoded in the overdensity 
parameter Abn 98 ( 2 )- We confirmed that there are variations due to 
it, but these are not sufficiently large to explain all of the scatter. 

d) initial mass function: As for different cosmologies and halo 
mass definitions, the assumption of different IMFs can be com¬ 
pensated for when calibrating the model: the observational data set 
used for the calibration (and the initial mass function assumed in its 
preparation) will determine the values of the model parameters - 
whatever the assumption for the model IMF. However, the model’s 
stellar masses (and other quantities not studied here such as the 
recycled fraction of gas, the amount of energy available for super¬ 
novae, chemical enrichment, etc.) will certainly be affected by the 
IMF choice. We can confirm that this has an effect but is not the 
primary source of the variations between models. 

e) tuning: This is by far the most decisive factor for the scatter 
(Henriques et al. 2009; Mutch et al. 2013). It has been shown by, 
for instance, Henriques et al. (2009) that DLB07 and an earlier 
version of the SAGE model could be brought into better agreement 
with the observational data by re-tuning their parameters optimally. 
Further, Lee et al. (2014) have shown that differences in merger 
trees could be overcome by re-tuning the model parameters. This 
directly applies to the comparison presented here: all models have 
been designed and tested using different simulations and merger 
tree (structures), but were not allowed to re-adjust their parameters 
for this initial project. We have seen that this has a very strong 
influence on the stellar mass function and is potentially the main 
source of the scatter seen in the plots throughout this paper. Future 
papers in this series will investigate the degree to which the models 
can be brought into agreement, and the extent to which they still 
differ, once retuned to the same set of observational constraints. 

We deliberately did not include any comparison with 
luminosity-based properties as their calculation involves another 
layer of complexity, e.g. stellar population synthesis, dust extinc¬ 
tion, etc. However - as stated in Appendix A - some models are 
using luminosity-related quantities to constrain their model param¬ 
eters. Also, one should not neglect the additional difficulties en¬ 
countered when moving from intrinsic galaxy properties such as 
mass to directly observable quantities such as luminosity or colour. 
Conversely, the stellar mass of a galaxy is not directly observable 
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and hence any derivation of it relies on modelling itself; therefore - 
as highlighted a couple of times before already - all observational 
data comes with its own error estimates that can be as large as 0.5 
dex at 2 = 0. All of this will certainly leave its imprint on the 
models presented here. 


8 CONCLUSIONS 

We conclude that applying galaxy formation models without due 
consideration to calibration with respect to cosmology, resolution, 
and - most importantly - merger tree prescription leads to scat¬ 
ter that could otherwise be avoided (see, e.g. Fontanot et al. 2009; 
Diaz-Gimenez & Mamon 2010; De Lucia et al. 2011; Contreras 
et al. 2013; Lu et al. 2014; Lee et al. 2014). But the need for re¬ 
calibration should not be viewed as a flaw of the models: it is a 
necessary step required to match the model to the particular obser¬ 
vational data sets that are chosen to underpin the model. The fact 
that a good match can be obtained is itself a non-trivial success of 
the model which indicates that the models capture the underlying 
key physical phenomena correctly. The (adjusted) parameters then 
place bounds upon the relevant physics, and the models can be used 
to test astrophysics outside that used in the calibration step. 

We close by mentioning again that this work only forms the 
initial step in a wider and long overdue programme designed to 
inter-compare current SAM and HOD models. The next stage is 
to calibrate all the models to a small, well specified, set of train¬ 
ing data, such as for instance the stellar mass function at 2 = 0 
and 2 = 2 before re-comparing the models on the other statistics 
shown in this work. This approach will likely significantly narrow 
the spread of the returned data on these physical quantities. It will 
also allow a more detailed comparison on such observationally in¬ 
teresting measures as hot and cold gas fractions, gas metallicity, 
galaxy sizes and morphologies, etc. The work for this has been 
started and will form the basis of a future workshop. 
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APPENDIX A: GALAXY FORMATION MODELS 
Al Galacticus (Benson) 

The Galacticus semi-analytic galaxy formation code (Benson 
2012) was designed to be highly modular. Every physical pro¬ 
cess is implemented through a simple and well-defined interface 


into which an alternative implementation of a calculation can eas¬ 
ily be added. Similarly, the physical description of galaxies is ex¬ 
tremely flexible. Each galaxy has a set of components (e.g. disk, 
dark matter halo, super-massive hlack hole, etc.) which can he cre¬ 
ated/destroyed as needed, each of which has a set of properties. 
The formation and evolution of galaxies is treated by simply defin¬ 
ing a set of differential equations for each galaxy. These are all 
simply fed in to an ordinary differential equation (ODE) solver 
which evolves them to a specified accuracy, removing any need 
for fixed timesteps. When running on merger trees derived from 
W-hody simulations, halo masses are interpolated linearly in time 
between available snapshots. In addition to this differential evo¬ 
lution galaxy components can define ‘interrupts’ so that the ODE 
solver stops, allowing the creation of new components (e.g. the first 
time gas cools and infalls it needs to create a disk component) or 
to handle discrete events (e.g. if a merger occurs the ODE solver is 
interrupted, the merger processed, and then the ODE solver starts 
up again). 

Cooling Cooling rates from the hot halo are computed using the 
traditional cooling radius approach (White & Frenk 1991), with a 
time available for cooling equal to the halo dynamical time, and 
assuming a /J-model profile with isothermal temperature profile (at 
the virial temperature). Metallicity dependent cooling curves are 
computed using CLOUDY (vl3.01, Ferland et al. 2013) assuming 
collisional ionization equilbrium. 

Star formation Star formation in disks is modeled using the pre¬ 
scription of Kmmholz et al. (2009), assuming that the cold gas of 
each galaxy is distributed with an exponential radial distribution. 
The scale length of this distribution is computed from the disks an¬ 
gular momentum by solving for the equilibrium radius within the 
gravitational potential of the disk-l-bulge-l-dark matter halo system 
(accounting for adiabatic contraction using the algorithm of Gnedin 
et al. 2004). 

Initial mass function A Chabrier (2003) IMF is used throughout. 

Metal treatment Metal enrichment is followed using the instanta¬ 
neous recycling approximation, with a recycled fraction of 0.46 and 
yield of 0.035. Metals are assumed to be fully mixed in all phases, 
and so trace all mass flows between phases. Metals affect the cool¬ 
ing rates from the hot halo, and also the star formation rate in disks 
(Krumholz et al. 2009). 

Supernova feedback and winds The wind mass loading factor, /3, is 
computed as /3 = (Vdisk/250km/s)“^'® where Vdisk is the circular 
velocity at the disk scale radius. Winds move cold gas from the disk 
back into the hot halo, where it remains in an outflowed phase for 
some fime before being reincorporated and possibly cooling once 
again. In the case of satellite galaxies, the ouflowing gas is added 
to the hot halo of the satellite’s host. 

Gas ejection & reincorporation Gas removed from galaxies by 
winds is retained in an outflowed reservoir. This reservoir gradu¬ 
ally leaks mass back into the hot halo on a timescale of fdyn/5 
where tdyn is the dynamical time of the halo at the virial radius. 

Disk instability The Efstathiou et al. (1982) criterion is used to 
judge when disks are unstable, with a stability threshold that de¬ 
pends on the gas fraction in the disk (0.7 for pure gas disks, 1.1 
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for pure stellar, linearly interpolated in between). When a disk is 
unstable, it begins to transfer stars and gas from the disk to the 
bulge on a timescale that equals the disk dynamical time for a max¬ 
imally unstable disk, and increases to infinite timescale as the disk 
approaches the stability threshold. 

Starburst There is no special “starburst” mode in Galacti- 
CUS. Instead, gas in the spheroid forms stars at a rate M* = 
0.04Mgaa/fdyn(V’/200km/s)“^, where fdyn is the dynamical 
time of the spheroid at its half mass radius, and V its circular ve¬ 
locity at the same radius. Starburst-level star formation rates are 
reached if enough gas is deposited into the spheroid, such as hap¬ 
pens after a merger. 

AGN feedback The mass and spin of black holes are followed in 
detail, assuming black holes accrete from both the hot halo and 
spheroid gas. When black holes are accreting from an advection 
dominated accretion flow, we compute the power of the jet pro¬ 
duced by the black hole using the method of Benson &. Babul 
(2009). This jet power is used to offset the cooling luminosity in 
the hot halo (if, and only if, the hot halo is in a hydrostatic phase), 
thereby reducing the cooling rate onto the galaxy. Additionally, a 
radiatively-driven wind is launched from the spheroid by the black 
hole, assuming that a fraction 0.0024 of its radiative output couples 
efficiently to the outflow. 

Merger treatment A merger between two galaxies is deemed to be 
“major” if their mass ratio exceeds 1:4. In major mergers, the stars 
and gas of the two merging galaxies are rearranged into a spheroidal 
remnant. In other, minor mergers, the merging galaxy is added to 
the spheroid of the galaxy that it merges with, while the disk of that 
galaxy is left unaffected. 

Substructures Substructures are traced using the subhalo informa¬ 
tion from the W-body simulation. 

Orphans When a subhalo can no longer be found in the W-body 
merger trees, a “subresolution merging time” is computed for the 
subhalo (based on its last known orbital properties and the algo¬ 
rithm of Boylan-Kolchin et al. (2008) algorithm. The associated 
galaxy is then an orphan, which continues to evolve as normal (al¬ 
though we have no detailed knowledge of its position within its host 
halo) until the subresolution merging time has passed, as which 
point it is assumed to merge with the central galaxy of its host halo. 

Calibration method The parameters of galaxy formation physics in 
Galacticus have been chosen by manually searching parameter 
space and seeking models which provide a reasonable match to a 
variety of observational data, including the z = 0 stellar mass func¬ 
tion of galaxies (Li & White 2009), 2 = 0 K and bj-band luminos¬ 
ity functions (Cole et al. 2001; Norberg et al. 2002), the local Tully- 
Fisher relation (Pizagno et al. 2007), the color-magnitude distribu¬ 
tion of galaxies in the local Universe (Weinmann et al. 2006), the 
distribution of disk sizes at z = 0 (de long & Lacey 2000), the 
black hole mass-bulge mass relation (Haring & Rix 2004), and the 
star formation history of the Universe (Hopkins 2004). GALACTI¬ 
CUS has also been calibrated to the local stellar mass function using 
MCMC techniques (Benson 2014), but so far only for a simplified 
implementation of galaxy formation physics. As such, we do not 
use this MCMC-calibrated version of GALACTICUS here. 


Model origin The parameters used were calibrated using Monte 
Carlo trees built using the algorithm of Parkinson et al. (2008). 

Modifications to the supplied data When importing merger trees 
for this project, GALACTICUS aims to make minimal changes to the 
tree structure. Two small modifications are required to ensure con¬ 
sistency of the merger trees. First, where halo A is indicated as be¬ 
ing the host of halo B, but A is not present in the merger trees (i.e. 
is not listed as the progenitor of any other halo), then A is assumed 
to be a progenitor of the same halo of which B is a progenitor. Sec¬ 
ond, if two haloes are mutual hosts (i.e. A is the host of B, while 
B is the host of A), GALACTICUS resolves this inconsistency by 
reassigning the more massive of the two haloes to be unhosted (i.e. 
to no longer be classified as a subhalo). No further modifications 
are made. In particular, this means that subhalo-subhalo mergers 
are allowed, subhaloes are allowed to become non-subhaloes later 
in their evolution, halo masses are permitted to decrease with time 
if indicated by the A^-body simulation, and subhaloes are permitted 
to jump between branches of merger trees (and between separate 
merger trees) if indicated. 

Halo finder properties used The standard incarnation of 
Galacticus uses MbnU but runs have been performed for 
all of the five supplied mass definitions. GALACTICUS further 
uses the following information from the provided halo catalogues: 
haloid, hosthaloid, number of particles, mass, radius, concentra¬ 
tion, spin parameter, angular momentum, position, and velocity; if 
any of these values is not supplied (e.g. spin parameter), a random 
value is drawn from a distribution as measured by cosmological 
simulations. Further, the peak value of the circular velocity curve 
and the velocity dispersion are carried through the code, but not 
used for any calculation. 


A2 GalICS- 2.0 (Cattaneo, Blaizot, Devriendt & Mamon) 

GalICS- 2.0 is not a simple development of GalICS (Hatton 
et al. 2003; Cattaneo et al. 2006), but a totally new and different 
code. Its main characteristics are presented here for the first time 
and hence described in more detail for than the other models. 

Gas accretion A baryonic mass Mb = /bMh is assigned to a halo 
of mass Mh. Let Treio be the temperature at which the intergalac- 
tic medium is reionised; /b = /b(Mh, z) is a function such that 
/b ~ 0 at Tvir < Treio and /b ~ Hb/GM at Tvlr > Feio- Its 
precise form is irrelevant for this article because the A-body sim¬ 
ulation used for this comparison has such poor resolution that it 
can resolve only haloes with Tvir S> Tieio- The baryonic mass that 
accretes onto a halo between two timesteps is the maximum be¬ 
tween AMb and zero. A fraction /hot = /hot(Mh, z) of this gas 
is shock heated and added to the hot halo. The rest is put into cold 
streams, which accrete onto the disk of the central galaxy on a dy¬ 
namical timescale. Cosmological hydrodynamic simulations show 
that /hot ~ 0 at Mh ^ 3 x 10^° Mq and that /hot ~ 1 at Mh > 
3x 10^^ Mq (Ocvirket al. 2008; Nelson et al. 2013). The halo mass 
M/i at which /hot = 0.5 increases with redshift (at least at 2 > 2). 
In this paper, we assume that fh is a linear ramp between the halo 
mass Mshock min at which some gas begins to be shock heated and 
the halo mass Mshock max above which the infalling gas is entirely 
shock heated. We assume that Mshock max = Mshock + criz — zf) 
for z > Zc and that Mshock max = Mshock for 2 ^ 2c where 
2c, Mshock and Mshock min are free parameters of the model to be 
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determined by fitting observational data. On a theoretical ground, 
A^shock is the critical mass halo mass at which the cooling time 
equals the gravitational compression time and Zc is the critical red- 
shift at which Mohock equals the non-linear mass (Dekel & Bim- 
boim 2006). 

Cooling Following Cattaneo et al. (2006) and Dekel et al. (2009), 
we assume that galaxies are built through the accretion of cold gas 
and that hot gas never cools. This is an extreme assumption, but 
is in good agreement with the galaxy colour-magnitude distribu¬ 
tion, while, if we let the gas cool, the predictions of the model 
are in complete disagreement with the observations (Cattaneo et al. 
2006). Introducing cooling makes sense only if one has a physi¬ 
cal model of how AGN feedback suppresses it. Attempts in this 
direction have been made, starting from (Croton et al. 2006), but 
the physics are uncertain. Hence it was considered that not much 
would be gained from implementing them in GalICS-2.0. 

Star formation Following (Bigiel et al. 2008), we assume a con¬ 
stant star formation timescale of = 2.5 x 10® yr for disks and 
f*b = 2.5 X 10® yr for bulges (merger-driven starbursts) for all 
gas. Mi, — Mgas/f*, where Mgas is the gas mass in the compo¬ 
nent. Star formation is suppressed when Sgaa < 1 Mq pc“®. Disks 
with Egaa > 20 Mq pc“® are assumed to be in a starburst mode 
and are assigned the same star formation timescale as bulges. Not 
only are these assumptions observationally motivated, but also they 
have the practical advantage that star formation rates are affected 
only mildly by errors in the modelling of disk sizes. 

Initial mass function As in Cattaneo et al. (2006), we assume a 
Kennicutt (1983) stellar initial mass function. The IMF that our 
model uses is determined by our choice of input stellar evolution 
tables (see Metal treatment below). 

Metal treatment Stars are created with metallicity that equals that 
of the gas from which they form. Stellar evolution is computed fol¬ 
lowing Pipino et al. (2009), who tabulated the mass loss rate and 
the metal yield of a stellar population as a function of its age and 
metallicity. We can select different stellar evolution models/IMFs 
by replacing the input file with the stellar evolution tables. 

Supernova feedback and winds Matter ejected from stars is mixed 
into the interstellar medium. The mass outflow rate through stel¬ 
lar feedback is assumed to be given by Mw = (Mh/MsN)^Af*, 
where M* is the star formation rate, while Msn (the halo mass at 
which outflow rate equals star formation rate) and j3 < Q are free 
parameters of the model. We consider a scaling with Mh rather 
than with Vc because the latter would lead to weaker feedback at 
high redshift for a same halo mass. 

Gas ejection & reincorporation As our model does not include 
cooling, it is irrelevant whether the gas that is blown out of the 
galaxy escapes from the halo or whether it is mixed to the hot com¬ 
ponent, as we assume currently. 

Disk structure Gas that accretes onto a galaxy contributes to the 
disk’s angular momentum Jd with specific angular momentum that 
equals that of the halo at the time of accretion. The disk radius is 
computed by assuming that the disk is exponential and by solving 


the equation 

A/f ^ 1 

|Jd| = ^ / e~M[vl{r)+ vl{r)+ vl{r)]^r^ dr, (Al) 

''"d Jo 

where the disk term Ud, the bulge term v^,, and the halo term Uh are 
the three terms that contribute in quadrature to disk’s rotation curve. 
The halo term Vh is computed by assuming that the halo follows an 
NFW profile modified by adiabatic contraction (Blumenthal et al. 
1986). 

Disk instability Disk instabilities are not enabled yet because we 
first want the test the properties of disks and bulges in a scenario in 
which mergers are the only mechanism for the formation of bulges. 

Starburst Gas is starbursting after a major merger or when 
Egaa > 100 Mq pc“®. In these cases, we lower the star formation 
timescale from = 2.5 x 10® yr to f*b = 2.5 x 10® yr. 

AGN feedback Two types of AGN feedback have been suggested 
in the literature: one is related powerful AGN, the other is mainly 
linked to Fanaroff-Riley type I radio sources. The need for the for¬ 
mer is unclear while the second is essential to explain the absence 
of cooling flows in massive systems (Cattaneo et al. 2009). Our 
code does not contain any explicit model of AGN feedback. How¬ 
ever, the assumption that the hot gas never cools is an implicit 
model for AGN feedback. It corresponds to the assumption that 
Tjet = T/x, i.e., that that jets self-regulate so that the power they 
damp into the hot gas matches exactly that which is lost to X-rays 
(see, e.g., Cattaneo & Teyssier 2007). Observationally, these quan¬ 
tities are equal to ~ 10% (Cattaneo et al. 2009, and references 
therein). Hence, this is a reasonable first approximation. 

Merger treatment In minor mergers (mass ratio < 1 : 3), the disk 
and the bulge of the smaller galaxy are added to the disk and the 
bulge of the larger galaxy, respectively. In major mergers, the galax¬ 
ies are scrambled into one large bulge and a fraction e, of their gas 
content feeds the growth of a central supermassive black hole (Cat¬ 
taneo et al. 2009). The size of the merger remnant is computed 
by applying an energy conservation argument that has been tested 
in hydrodynamic simulations (Covington et al. 2011; Oser et al. 
2012). Our calculation assumes that merging pairs start with zero 
energy of interaction at infinity in agreement with what we see in 
cosmological hydrodynamic simulations and with the constraints 
from the mass-size relation (Shankar et al. 2013, 2014). 

Substructures Substructures are traced from the A^-body simula¬ 
tion. No gas accretion is allowed onto them but cold gas already 
accreted keeps streaming onto satellite galaxies as long as the host 
halo mass is Mh < Mshockmax- 

Orphans Whenever the code encounters a halo with more than one 
progenitor, it computes the dynamical friction time for all progeni¬ 
tors bar the most massive one. The dynamical friction time is com¬ 
puted with Jiang et al. (2008)’s formula as in Cattaneo et al. (2011). 
This formula is a modification of Chandrasekhar (1943)’s (see Bin- 
ney & Tremaine 2008) that includes the effects of orbital eccen¬ 
tricity and that has been calibrated on the results of cosmological 
simulations. Our calculation of the dynamical friction time con¬ 
tains a fudge factor Cdf that is a free parameter of the model. A 
halo/subhalo merges with its descendent according the merger tree 
only after this time has elapsed. The halo catalogues are completed 
with the creation of a ghost halo at all timesteps between the one 


(gl 2010 RAS, MNRAS 000, 1-34 


nlFTy Galaxies 23 


when the halo/subhalo was last detected and the elapsing of the dy¬ 
namical friction time. If Cdf = 0, the dynamical friction time is set 
equal to zero. Hence halo/subhalo mergers result into immediate 
galaxy mergers. 

Parameters and calibration method For this project, the code has 
been calibrated manually by requiring that it fits the evolution of 
the galaxy stellar mass function at 0 < 2 : < 2.5 (Baldry et al. 
2012; Ilbert et al. 2013; Bernard! et al. 2013). The parameters that 
have been tuned and that are relevant for this comparison are: 

• The mass Mshockmin above which shock heating begins and 
the three parameters Mshock, a and Zc that determine the critical 
mass above which shock heating is complete. 

• The star formation timescales for disks (f*d) and bulges 
(merger-driven starbusts; f*d), and the gas surface density thresh¬ 
old for star formation Emin. 

• The mass Msn at which outflows rate equal star formation 
rate and the exponent jd of the scaling of mass-loading factor with 
halo mass. 

• The dynamical friction parameter Cdf. 

• The critical mass ratio /r that separates minor and major merg¬ 
ers. 

Here we have used Afghockmin — 10 ATq, ATghock — 

M 0 , a = 0.3, Zc = 1.5 and Msn = 4 x lO’^^ Mq (for 
h = 0.678). The best fit to the galaxy stellar mass function of 
Baldry et al. (2012) is found for /3 = —2 and Cdf > 0 (our best 
fit is for Cdf ^ 1; hence we apply Occam’s razor and prefer a 
model with no ghost/orphans). The star formation timescales for 
disks (f*d = 2.5 X 10® yr) and merger-driven starbursts (f*d = 
2.5 X 10* yr) and the gas surface density threshold for star forma¬ 
tion Emin = 1 Mq pc“® Were not tuned but they were set to the 
values found observationally by Bigiel et al. (2008). The critical 
mass ratio of /i = 1/3 that separates minor and major mergers was 
not tuned either. 

Model origin The model was conceived to be run on merger trees 
from A^-body simulations. 

Modifications to the supplied data None besides the creation of 
ghost haloes if ear > 0. 

Halo finder properties used They halo properties that enter the 
GalICS- 2.0 SAM are the halo mass Mhaio, its radius i?haio and 
the halo angular momentum Jhaio (the virial velocity is computed 
with the formula n/ = GMhaio/7?haio). Positions and velocities 
are used to compute dynamical friction and tidal stripping. Halo 
concentration is needed to compute dynamical friction, tidal strip¬ 
ping, and the radii of disks and bulges, but it is not provided in the 
SUBFIND catalogues. Hence, it has been computed with the fit¬ 
ting formulae of Munoz-Cuartas et al. (2011). GalICS- 2.0 flags 
certain haloes as bad. Bad haloes have positive total energy or they 
are tidal features that broke off a good halo, in which case we as¬ 
sume that no galaxy is formed in them. Normally, there is no gas 
accretion onto haloes that are flagged as bad. We could not im¬ 
pose this condition here because halo kinetic and potential ener¬ 
gies were not provided but this is unlikely to affect our results be¬ 
cause the halo/subhalo finder SUBFIND automatically removes 
unbound particles. 


A3 Morgana (Monaco & Fontanot) 

In this paper we use the version of MORGANA that has been pre¬ 
sented in Lo Faro et al. (2009). Its main properties are detailed be¬ 
low. 

Cooling Described in detail in Viola et al. (2008), the model as¬ 
sumes that the gas is in hydrostatic equilibrium within an NFW 
halo and polytropic. Its thermal state responds to the injection of 
energy. The cooling radius is treated as a dynamical variable. Disk 
sizes are computed using the Mo et al. (1998) model, taking into 
account the presence of a bulge. 

Star formation Star formation is treated following the results of the 
model by Monaco (2004). 

Initial mass function A Chabrier (2003) IMF is assumed. 

Metal treatment An Instantaneous Recycling Approximation is as¬ 
sumed, only global metallicities are followed. 

Supernova feedback and winds Feedback and winds are treated 
following the results of the model by Monaco (2004). The ejec¬ 
tion rate of gas into the halo is always equal to the star formation 
rate. 

Gas ejection & reincorporation The halo acts as a buffer for feed¬ 
back. Galaxies inject (hot or cold) mass and (thermal or kinetic) 
energy in the halo, when the typical specific energy is larger than 
the escape velocity of the halo the gas is ejected into the IGM. Half 
of this ejected gas is re-accreted whenever the circular velocity of 
the halo grows larger than the velocity at ejection. 

Disk instability The Efstathiou et al. (1982) criterion is used. 

Starburst The size of star-forming gas in bulges is estimated by 
assuming that its velocity dispersion is determined by turbulence 
injected by SN feedback. SFR is computed using a Kennicutt law, 
the high gas surface densities guarantees short gas consumption 
timescales. 

AGN feedback Star formation in the bulge is responsible for the 
loss of angular momentum of a fraction of gas; the reservoir of such 
low-angular momenum gas is accreted onto the BH on a viscous 
time-scale. The accretion is Eddington-limited. AGN feedback can 
be in the quasar mode or in the ‘jet’ (radio) mode if the accretion 
rate in units of the Eddington one is higher or lower than 0.01. 

Merger treatment Major mergers transform the whole resulting 
galaxy into a bulge, a minor merging satellite is put in the galaxy 
bulge. A fraction of 80% of stars in minor merging satellites is po¬ 
sitioned into the stellar halo component. 

Substructures The W-body simulation subhalo information is not 
used to follow the orbital evolution of substructures; they are rather 
tracked analytically. In detail, whenever a DM halo merges with a 
larger structure, the orbital decay of its galaxies is computed using 
an updated version of the fitting formulae provided by Taffoni et al. 
(2003), calibrated on the results of numerical simulations, which 
account for dynamical friction, tidal stripping and tidal disruption. 
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Orphans All satellite galaxies are effectively treated as ‘orphans’, 
i.e. the subhalo information is not explicitly used in modeling their 
evolution. 

Calibration method Parameters were manually to fit the stellar 
mass function of galaxies at 2 = 0 and the evolution of the SFR 
density; see Lo Faro et al. (2009) for more details. 

Model origin MORGANA has been designed to work with the ide¬ 
alized merger trees obtained with the PiNOCCHIO code (Monaco 
et al. 2002), so the application to numerical merger trees requires a 
significant amount of cleaning of the trees. 

Modifications to the supplied data It is assumed that the merger 
trees of haloes are mirrored by the evolution of their main haloes, 
so that main haloes either have no progenitors or have at least one 
main halo progenitor. When this is not true, the merger tree is mod¬ 
ified to adapt to this requirement. This includes the reabsorption 
of any substructure that descends into a main halo (a ‘backsplash’ 
halo) or the exchange of the role of main halo between two haloes. 

Halo finder properties used As a halo mass, our adopted default 
choice is the FOF mass Mfof, but we run the model on all five mass 
definitions. However, the halo mass used to obtain the budget of 
baryons available to the halo is assumed never to decrease, so it is 
computed as the maximum of the total mass that the halo and all 
its progenitors have got in the past. In addition to this MORGANA 
made use of the provided haloid, hosthaloid, number of substruc¬ 
tures, number of particles, position, and velocity. 

A4 Sag (Cora, Vega-Martmez, Gargiulo & Padilla) 

While the SAG model originates from a version of the Munich code 
(Springel et al. 2001b), it has seen substantial development and im¬ 
provement; we are going to explain its features (and derivations 
from the Munich and other models) here. 

Cooling A gaseous disk with an exponential density profile is 
formed from gas inflow generated as the result of the radiative gas 
cooling suffered by the hot gas in the halo (Springel et al. 2001). 
The metal-dependent cooling function is estimated by consider¬ 
ing the radiated power per chemical element given by Foster et al. 
(2012) multiplied by the chemical abundances of each element in 
the hot gas, thus being completely consistent with the metallicity 
of this baryonic component. 

Star formation When the mass of the disk cold gas exceeds a crit¬ 
ical limit, an event of quiescent star formation takes place, as in 
Croton et al. (2006). 

Initial mass function SAG assumes a Salpeter IMF. 

Metal treatment SAG includes a detailed chemical implementation 
(Cora 2006), which estimates the amount of metals contributed by 
stars in different mass ranges. Metals are recycled back to the cold 
gas taking into account stellar lifetimes (Padovani & Matteucci 
1993). The code considers yields from low- and intermediate-mass 
stars (Karakas 2010), mass loss of pre-supernova stars (Hirschi 
et al. 2005), and core collapse supernovae (SNe CC, Kobayashi 
et al. 2006). Ejecta from supemovae type la (SNe la) are also in¬ 
cluded (Iwamoto et al. 1999); SNe la rates are estimated using the 
single degenerate model (Lia et al. 2002). 


Supernova feedback and winds Feedback from SNe CC is mod¬ 
eled following De Lucia et al. (2004a). The amount of reheated 
gas is proportional to the energy released by SNe CC and inversely 
proportional to the square of the halo virial velocity. Cora (2006) 
adapted this prescription according to the chemical model imple¬ 
mented, such that the energy contribution of SNe CC occurs at 
the moment of their explosions, for which the lifetimes of massive 
stars are considered. SNe feedback produces outflows of material 
that transfer the reheated cold gas with its metal content to the hot 
gas phase; this chemical enrichment has a strong influence on the 
amount of hot gas that can cool, since the cooling rate depends on 
the hot gas metallicity. 

Gas ejection & reincorporation SAG assumes a ‘retention’ scheme 
in which the cold gas reheated by SNe feedback is transfered to the 
hot gas phase, being available for further gas cooling that takes 
place only in central galaxies (‘strangulation’ scheme). 

Disk instability When a galactic disk is sufficiently massive that its 
self-gravity is dominant, it becomes unstable to small perturbations 
by satellite galaxies. The circular velocity of the disk involved in 
the stability criterion of Cole et al. (2000) is approximated by the 
velocity calculated at three times the disk scale-length, where the 
rotation curve flattens (see Tecce et al. 2010, for details concerning 
disk features). We model the influence of a perturber by comput¬ 
ing the mean separation between galaxies in a group; the instability 
is triggered when this separation is smaller than a certain factor (a 
free parameter of the model) of the disk scale radius of the per¬ 
turbed galaxy. When an unstable disk is perturbed, existing stars 
are transferred to the bulge component along with the cold gas that 
is consumed in a starburst. 

Starburst Starbursts occur in both mergers and triggered disk in¬ 
stabilities and are the only channel for bulge formation. During 
a starburst episode, we consider that the cold gas available in the 
bulge that has been transfered from the disk, referred to as bulge 
cold gas, is gradually consumed. The period for star formation is 
chosen to be the dynamical time-scale of the disk. However, as the 
starburst progresses, effects of supemovae feedback, recycling of 
gas from dying stars and black hole growth modify the reservoir of 
cold gas of both disk and bulge, thus also changing the time-scale 
of the starburst (Gargiulo et al. 2014). 

AGNfeedback SAG includes radio-mode AGN feedback following 
Croton et al. (2006) as described in Lagos et al. (2008), which re¬ 
duces the amount of gas that can cool thus providing a mechanism 
for regulating star formation in massive galaxies. AGNs are pro¬ 
duced from the growth of central black holes, for which two possi¬ 
ble mechanisms are considered: i) infall of gas towards the galactic 
centre, induced by merger events or disk instabilities, and ii) the 
accretion of gas during the cooling process. The current version of 
the code (Gargiulo et al. 2014) considers that mass accretion rates 
in the latter case depends on the square of the virial velocity (being 
consistent with a Bondi-type accretion, Bondi 1952), instead of on 
the cube of the velocity as in Lagos et al. (2008), in order to prevent 
super massive black holes at the centre of cluster-dominant galax¬ 
ies to grow unrealistically large, at the expense of the intracluster 
medium. 

Merger treatment The galaxy inhabiting the subhalo is assumed to 
merge with the central galaxy of their host subhalo on a dynami¬ 
cal friction time-scale (Binney & Tremaine 1987). The merger can 
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be major or minor depending on the the baryonic mass ratio be¬ 
tween the satellite galaxy and the central galaxy. If this mass ratio 
is larger than 0.3, then the merger is considered a major one. In this 
case, stars and cold gas in the disk of the remnant galaxy are trans¬ 
ferred to the bulge, with the latter being consumed in a starburst. 
The presence of a starburst in a minor merger will depend on the 
fraction of cold gas present in the disk of the central galaxy, as im¬ 
plemented by Lagos et al. (2008). In minor mergers, only the stars 
of the merging satellite are transferred to the bulge component of 
the central galaxy. 

Substructures Substructures are followed from the A^-body simu¬ 
lation. The supplied data has not been modified. Those branches of 
merger trees that start with subhaloes give place to spurious galax¬ 
ies with neither cold gas nor stars since gas cooling does not take 
place in galaxies residing in subhaloes. Those merger trees are ig¬ 
nored and we discard those galaxies in the output of SAG. 

Orphans Orphan galaxies are created when their subhaloes are no 
longer identified after merging with a larger one because they lose 
mass as a result of tides. The position and velocity of orphan galax¬ 
ies are estimated assuming a circular orbit with a velocity given by 
the virial velocity of the parent subhalo and a decaying radial dis¬ 
tance determined by dynamical friction (Binney & Tremaine 1987). 

Calibration method Calibrations of SAG are performed using the 
‘Particle Swarm Optimization’ technique, which yields best-fitting 
values for the free parameters of the model allowing it to achieve 
good agreement with specific observational data (Ruiz et al. 2015). 
The free parameters that have been tuned are those related with star 
formation efficiency, the SNe feedback efficiencies that control the 
amount of disk cold gas and bulge cold gas reheated by the energy 
generated by SNe formed in the disk and in the bulge, respectively, 
parameters involved in the AGN feedback, that is, the fraction of 
cold gas accreted onto the central supermassive black hole (SMBH) 
and the efficiency of accretion of hot gas onto the SMBH, and fi¬ 
nally, the factor involved in the distance scale of perturbation to 
trigger disk instability. 

We calibrate the free parameters of the SAG model consid¬ 
ering a set of observational constraints that involve the 2 = 0 lu¬ 
minosity function in the r-band (r-band LF), the relation between 
the mass of the central SMBH and the bulge mass (BHB relation), 
the redshift evolution of SNe la and SN CC rates, and the [a/Fe]- 
stellar mass relation of elliptical galaxies. The first two constraints 
help to tune the free parameters associated to the star formation ef¬ 
ficiency and the SNe and AGN feedback. The third one allows to 
fix the fraction of binary stars that explode as SNe la, and therefore, 
the amount of iron recycled into the interstellar medium. The last 
one establishes more restriction to the efficiency of SNe feedback 
arising from stars formed in the bulge. The observational data used 
are the r-band LF of Blanton et al. (2005a), the BHB relation given 
by Haring & Rix (2004) and Sani et al. (2010), the compilation of 
rates for both SNe la and SNe CC given by Melinder et al. (2012), 
and the [a/Fe] ratio of elliptical galaxies presented in Thomas et al. 
(2010) and Arrigoni et al. (2010). 

Model origin While the SAG model originates from a version of 
the Munich code (Springel et al. 2001), based on A^-body simula¬ 
tions, it has seen substantial development and has been improved 
with a detailed chemical implementation (Cora 2006), the inclu¬ 
sion of AGN feedback and disk instabilities (Lagos et al. 2008), a 
detailed estimation of disk scale-lengths Tecce et al. (2010), and 


a gradual star formation during starbursts (Gargiulo et al. 2014), 
among other aspects that are not taken into account in the version 
used for the current comparison, like the effects of accretion with 
misaligned angular momenta on the properties of galactic disks 
(Padilla et al. 2014), a star formation dependent top-heavy inte¬ 
grated galactic IMF (Gargiulo et al. 2014), estimation of nebular 
emission of star-forming galaxies (Orsi et al. 2014), and environ¬ 
mental effects such as tidal stripping and ram pressure stripping 
(Cora et al., in prep.). 

Modifications to the supplied data As mentioned above, branches 
of the merger tree that start with subhaloes are ignored. However, 
they contribute to generate the galaxy population if gradual removal 
of hot gas is allowed since, in that case, satellites receive cooling 
flows (which, however, is not the case in the version of the model 
considered here). No other modifications have been made. 

Halo finder properties used SAG uses the bound mass Mbnd to 
construct its galaxies. Further properties entering SAG are the 
haloid, hosthaloid, number of substructures, radius, position, ve¬ 
locity, and spin parameter. 

AS SantaCruz (Somerville & Hirschmann) 

The SantaCruz SAM was first presented in Somerville & Pri- 
mack (1999) and significantly updated in Somerville et al. (SOS 
2008), Somerville et al. (2012), Hirschmann et al. (2012, H12), 
and recently in Porter et al. (2014, P14). The SantaCruz SAM 
includes the following physical processes: (1) shock heating and 
radiative cooling of gas, (2) conversion of cold gas into stars via an 
empirical ‘Kennicutt-Schmidt’ relation, (3) starbursts, black hole 
feeding, and morphological transformation due to galaxy mergers 
4) bulge and black hole growth via ‘disk instabilities’ 5) metal en¬ 
richment of the interstellar and intracluster media by supemovae 
using the instantaneous recycling approximation 6) galactic out¬ 
flows driven by stars and supemovae (7) galactic outflows driven 
by ‘quasar mode’ black hole activity, and heating of the hot intra¬ 
cluster and intragroup medium by ‘radio mode’ black hole activity. 
The code used here adopts the modifications relative to the 808 
model suggested by Hirschmann et al. (2012) in order to match the 
observed luminosity function of AGN: massive seed black holes, 
black hole feeding via disk instabilities, and suppressed black hole 
feeding in mergers at low redshift (see H12 for details). These mod¬ 
ifications have almost no impact on galaxy properties but do affect 
black hole growth. 

Cooling The cooling model is based on the spherically symmetric 
cooling flow model originally presented in White & Frenk (1991), 
and is described in detail in 808 and P14. We assume a singular 
isothermal profile for the hot gas density distribution and adopt 
the metallicity dependent cooling functions of Sutherland & Do- 
pita (1993). If the cooling radius is larger than the virial radius, gas 
is accreted on a halo dynamical time. Cooling gas is only accreted 
onto ‘central’ galaxies. 

Star formation As gas cools, we assume that it settles into a 
rotationally-supported exponential disk. We use a prescription 
based on the work of Mo et al. (1998) and Somerville et al. (2008) 
to compute the radial size of the gas disk. 

We allow for two modes of star formation: ‘disk mode’ star 
formation, which occurs in disks at all times as long as cold gas 
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above a critical surface density is present, and ‘burst mode’ star for¬ 
mation, which occurs after two galaxies merge. In the ‘disk mode’ 
the star formation rate density is dependent on the surface density 
of cold gas in the disk, following the empirical Schmidt-Kennicutt 
relation (Kennicutt 1998). Only gas that is above a critical surface 
density threshold is allowed to form stars. See SOS and P14 for fur¬ 
ther details. 

Initial mass function We adopt a Chabrier (2003) initial mass func¬ 
tion. 


Metal treatment We model chemical enrichment using the instan¬ 
taneous recycling approximation. The chemical yield is treated as 
a free parameter. 


Supernova feedback and winds Massive stars and supemovae are 
assumed to produce winds that drive cold gas back into the ICM 
and IGM, heating the gas in the process. The mass outflow rate is 
proportional to the star formation rate. 


rhrh = esN 
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where csn and arh are free parameters and Vdisk is the circular 
velocity of the disk. 


Gas ejection & reincorporation The proportion of the gas that is 
ejected from the halo entirely is a decreasing function of the halo’s 
virial velocity. This gas can then fall back into the hot halo, at a 
re-infall rate that is proportional to the mass of the ejected gas and 
inversely proportional to the dynamical time of the halo (see 808 
and P14 for details). When gas is ejected due to supernova feed¬ 
back, these winds are assumed to have a metallicity Zcoid. Ejected 
metals are assumed to “re-infall” back into the hot halo on the same 
timescale as the gas. 


Disk instability We include bulge formation and black hole growth 
due to disk instabilities as in the “stars” model of P14 (see also 
H12). Following Efstathiou et al. (1982) we define the stability pa¬ 
rameter 


Cdisk — 
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where IGiax is the maximum circular velocity of the halo, ruisk is 
the scale length of the disk, and Mdiak is the mass of stars in the 
disk. Disks are deemed unstable if Cdisk < Ccrit, where Ccrit is a free 
parameter. In every timestep in which the disk becomes unstable, 
we move just enough stars from the disk to the bulge component to 
restore stability. 


Starburst Following a merger with mass ratio p > /dcrit ~ 0.1, 
we trigger a burst of star formation. The burst mode star formation 
is added onto the ‘disk’-mode star formation described above. The 
efficiency and timescale of the burst mode depends on the merger 
mass ratio, the gas fraction of the progenitors, and the circular ve¬ 
locity of the progenitors, as described in 808. These scalings were 
derived from hydrodynamic simulations of binary galaxy-galaxy 
mergers (see 808 for details). 


are assumed to merge fairly quickly, and the resulting hole grows at 
its Eddington rate until the energy being deposited into the I8M in 
the central region of the galaxy is sufficient to significantly offset 
and eventually halt accretion via a pressure-driven outflow. This 
results in self-regulated accretion that leaves behind black holes 
that naturally obey the observed correlation between BH mass and 
spheroid mass or velocity dispersion (see 808 and HI2 for more 
details). Large-scale winds associated with this rapid BH growth 
can also remove gas from the galaxy (see 808). 

A second mode of black hole growth, termed “radio mode”, 
is assumed to couple very efficiently with the hot halo gas, and 
to provide a heating term that can partially or completely offset 
cooling during the “hot flow” mode (we assume that the jets cannot 
couple efficiently to the cold, dense gas in the infall-limited or cold 
flow regime). 

Merger treatment Once haloes cross the virial radius of a parent 
halo, the 8anta Cruz SAM tracks the timescale for the orbits of 
the sub-haloes (satellites) to decay via dynamical friction using a 
refined version of the Chandrasekhar formula (see 808). Dark mat¬ 
ter is stripped from the sub-haloes on each orbit. If the sub-halo 
is stripped below a critical mass, it and the galaxy it contains are 
considered to be tidally disrupted and the stars in the galaxy are 
added to the ‘diffuse stellar halo’. Any cold gas is considered to be 
heated and added to the hot gas halo of the host. Galaxies that are 
not tidally destroyed before they reach the center are merged with 
the cenhal galaxy. Following a merger, a fraction /scat of the stars 
from the satellite are added to the diffuse interstellar halo. Depend¬ 
ing on the merger mass ratio and the gas fraction of the progenitors, 
some of the disk stars are moved to the bulge following a merger 
(see P14). 

Substructures Sub-structures are tracked analytically. One minor 
issue is that our SAM cannot track haloes that become sub-haloes 
and then travel outside the virial radius of the host (so-called ‘back- 
splash galaxies’, e.g. Gill et al. 2005)- we continue to treat these as 
sub-haloes. 

Orphans All substructures are effectively treated as ‘orphans’ in 
our models. 

Calibration method The SAM parameters used here are identical 
to those used in Porter et al. (2014). These parameters were cho¬ 
sen by tuning to observations of the galaxy stellar mass function, 
gas fraction and stellar metallicity as a function of stellar mass, 
and black hole mass versus bulge mass relationship, all at 2 ~ 0 
(see 808 and P14 for details). In addition, our “disk instability” pa¬ 
rameter is tuned to attempt to reproduce the morphological mix of 
nearby galaxies as a function of stellar mass (see P14). 

Model origin The SantaCruz SAM was originally developed 
based on EPS merger trees (Somerville & Primack 1999). Sub¬ 
sequently, the model has been implemented within merger trees 
extracted from the Bolshoi simulations (P14). Our model results 
are quite insensitive to whether we use EPS mergers trees or high- 
quality A^-body merger trees. 


AGN feedback In addition to triggering starbursts, mergers drive 
gas into galactic nuclei, fueling black hole growth. Every galaxy is 
born with a small “seed” black hole (typically ~ 10'* Mq in our 
standard models). Following a merger, any pre-existing black holes 


Modifications to the supplied data Currently the SantaCruz 
SAM does not make use of positional information from the A^-body 
simulation for “sub-haloes”, haloes that have become subsumed in 
other virialized haloes. As a result sub-haloes were stripped from 
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the merger trees provided for this project. No further modifications 
were made. 

Halo finder properties used The SantaCruz SAM uses the 
ATbn 98 ‘virial mass’ definition by default. The only halo finder 
properties that we use are the halo mass and redshift, and progenitor 
and descendant relationships. 


A6 ySAM (Lee & Yi) 

Here we briefly summarize the main features of ySAM. It has 
been developed to calculate galaxy properties on halo merger trees 
extracted from A^-body simulations. ySAM assumes that haloes 
newly identified in a volume have hot gas components propor¬ 
tional to their virial mass, following the universal baryonic fraction, 
flb/flm. The hot gas components are cooled (see below). Cold gas 
components form gas disks and stellar populations are newly bom 
in the disks by the simple law in Kauffmann et al. (1993). Grad¬ 
ual mass loss from stellar populations can enrich the metallicity 
of galaxies by recycling the ejecta in hot and cold gas reservoirs. 
ySAM also includes feedback processes (see below). In addition, 
environmental issues can affect the gas components of subhaloes: 
hot gas is stripped by tidal forces (see Kimm et al. 2011) and ram 
pressure (McCarthy et al. 2008; Font et al. 2008). Perhaps the most 
notable difference of yS AM is that it calculates stellar mass loss in 
all constituent stellar populations in each galaxy step by step. This 
often involves tracking tens of thousands of separate populations 
in each galaxy, which helps to trace the gas recycling more real¬ 
istically than in the case of instant recycling assumption. Further 
details can be found in Lee & Yi (2013). 

Cooling ySAM calculates gas cooling rates by adopting the pre¬ 
scription proposed by White & Frenk (1991). 

Star formation ySAM follows a simple law suggested by Kauff¬ 
mann et al. (1993) for quiescent star formation. 

Initial mass function Chabrier, Salpeter, and Scalo IMFs are avail¬ 
able in ySAM. Stellar populations formed in quiescent and bursty 
modes can have different IMFs. 

Metal treatment ySAM calculates stellar mass loss from every 
single stellar population at each time step. The amount of metals in 
the ejecta is computed based on a given IMF and chemical yields 
of Type la (Nomoto et al. 1984) and Type II SNe (Portinari et al. 
1998), and intermediate mass stars (Marigo 2001). The metals can 
be recycled by star formation or circulated between galaxies and 
environments via gas cooling or heating. 

Supernova feedback and winds ySAM follows the prescriptions 
described in Somerville et al. (2008) for SN feedback and winds. 

Gas ejection & reincorporation Gas components can be blown 
away by SN and QSO mode AGN feedback. Some of them can 
be re-accreted onto galaxies in the dynamical timescale of haloes. 

Disk instability In ySAM, disk instability can be estimated by us¬ 
ing a formula derived by Efstathiou et al. (1982). Due to uncertain¬ 
ties of disk instability (e.g. Athanassoula 2007), however, we turn 
it off in this study. 


Starburst The amount of stars bom in bursty mode is evaluated 
when galaxy mergers (M 2 /Mi > 0.1) take place. It is calculated 
following the prescriptions formulated by Somerville et al. (2008) 
from the numerical simulations performed by Cox et al. (2008). 

AGN feedback QSO and radio modes AGN feedback has been im¬ 
plemented into ySAM by following Croton et al. (2006). 

Merger treatment Satellite galaxies merge into their centrals when 
subhaloes harbouring them reach very central region of host haloes 
(< 0.1J(vir). If a subhalo is deprived of mass below baryonic mass, 
then the galaxy at the centre of the subhalo is considered to be dis- 
mpted and their stellar components become diffuse stellar compo¬ 
nents of its host halo. 

Substructures ySAM traces substructures following the results 
from the A-body simulation. 

Orphans If a substructure disappears before reaching the central 
region of its host halo, ySAM calculates its mass (Battin 1987) 
and orbit (Binney & Tremaine 2008) analytically until approaching 
the very central regions. This has a large impact on the lifetime of 
subhaloes and galaxy merging timescale (Yi et al. 2013). 

Calibration method yS AM has been manually calibrated to match 
galaxy mass functions (mainly that in Panter et al. 2008), BH-bulge 
relation (Haring & Rix 2004), global star formation density evolu¬ 
tion (Panter et al. 2008), and stellar-to-halo mass relation (Moster 
et al. 2010). 

Model origin ySAM has been developed to be run based on halo 
merger trees extracted from A-body simulations. 

Modifications to the supplied data We ignore halo merger trees 
that disappear as independent host haloes before merging into other 
haloes. We also remove merger trees identified as subhaloes at the 
beginning. If, however, haloes bom as subhaloes come out of their 
hosts and remain as host haloes by 2 = 0, then we do not discard 
them. 

Halo finder properties used ySAM adopts as its prime mass 
M 200 C and also uses the following information from the supplied 
halo catalogues: haloid, hosthaloid, radius, position, and velocity; 
the number of particles, peak value and position of the circular ro¬ 
tation curve, and spin parameter are also read, but not used. 

A7 Durham - Galform (Gonzalez-Perez, Bower & Font) 

Galform-gp14, Galform-kb06 & Galform-kf08 

For the study presented here we use the Gonzalez-Perez et al. 
(2014) (thereafter Galform-GP14) flavour of the Galform 
model (Cole et al. 2000), which exploits a Millenium Simulation- 
class A-body run performed with WMAP7 cosmology. We also 
compare the results from two variations of the Galform-GP14 
model, to which we refer to as GALFORM-KB06 and Galform- 
KF08. These have been generated by mnning a modified version 
of the Galform-GP14 that accounts for the main differences with 
respect to the Bower et al. (2006) and Font et al. (2008) models. 
These two models were developed using merger trees derived from 
the Millennium simulation (Springel et al. 2005), which assumes a 
cosmology close to that from WMAPl. 
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The model referred here as GALFORM-KB06 is the 
Galform-GP 14 model run assuming a single star formation law. 
The model referred here as Galform-KFOS is the Galform- 
GP14 model run assuming a single star formation law and including 
a gradual stripping of hot gas in satellite galaxies as opposed to the 
strangulation assumed by default. 

Cooling Cooling rates are estimated by defining a cooling radius, 
assuming that the shock-heated halo gas is in collisional ionization 
equilibrium. The gas density profile in the halo is kept fixed and is 
well fit by the /3-model (Cavaliere & Fusco-Femiano 1976) gener¬ 
ally used to model hot X-ray emitting gas in galaxy clusters. We 
use the cooling functions tabulated by Sutherland & Dopita (1993), 
which are a function of both the metallicity and temperature of the 
gas. 

Star formation The Galform-GP 14 model assumes two different 
star formation laws depending on the star formation being quies¬ 
cent, happening in disks, or happening in a burst. In the star burst 
mode, the star formation rate is assumed to be simply proportional 
to the mass of cold gas present in the galaxy and inversely propor¬ 
tional to star formation time-scale. This is what is also assumed 
for the quiescent star formation in both the GALFORM-KB06 and 
Galform-KF 08 model variants. The quiescent star formation in 
the Galform-GP 14 model is obtained in a self consistent calcu¬ 
lation in which HI and H 2 are tracked explicitly and the star forma¬ 
tion in disks is assumed to depend on the amount of molecular gas, 
H 2 , rather than on the total mass of cold gas (Lagos et al. 2011). 

IMF A Kennicutt IMF is assumed. 

Metal treatment The model uses the instantaneous recycling ap¬ 
proximation, with a recycled fraction of 0.44 and a metal yield of 
0 . 021 . 

Supernova feedback and winds The supernova feedback efficiency 
is quantified in the model in terms of the rate at which cold gas 
is reheated and thus ejected into the halo, Mreheated, per unit 
mass of stars formed, ft, which are computed as Mreheated = 



Gas ejection & reincorporation The gas affected by stellar feed¬ 
back is assumed to be heated to the virial temperature of the current 
halo and placed into a reservoir. The mass in this reservoir, Mrea, 
will return to the hot halo at a rate given by 1.26 x Mrea/fdyn, 
where fdyn is the dynamical time-scale of the halo. The gas that 
has been reincorporated into the halo can then cool back on to the 
galaxy disk. 

Disk instability If a disk is strongly self-gravitating it will be un¬ 
stable with respect to the formation of a bar. This will happen in 
those disks satisfying the Efstathiou et al. (1982) criterion, with a 
stability threshold of 0.8. 

Starburst During star bursts episodes, the star formation law as¬ 
sumed in the Galform-GP 14 model is different from the quies¬ 
cent case (not so for the other two model variants). The available 
cold gas is assumed to be consumed during a starburst event with a 
finite duration. 


AGN feedback The onset of the AGN supression of the cooling 
flow can only happen in the model in haloes undergoing quasi¬ 
hydrostatic cooling. This is assumed to happen for haloes hosting 
galaxies such that fcooi > fff/0.6, where tcooi is the cooling time 
of the gas and f// is the free-fall time for the gas to reach the centre 
of the halo. 

Merger treatment Mergers such that the ratio between masses ex¬ 
ceeds 0.1 will trigger a burst of star formation. Disks are trans¬ 
formed into spheroids when mergers happen. 

Substructures The A^-body simulation subhalo information is used 
to trace substmcture. 

Orphans When the subhalo hosting a satellite galaxy can no longer 
be followed with the N-body simulation information, the Lacey & 
Cole (1993) analytical expression is used to compute the merging 
timescale of this orphan galaxy. 

Calibration method The free parameters in this model where cho¬ 
sen manually such that the predicted luminosity functions in both 
bj and K-band at redshift 0 and the predicted evolution of the 
rest-frame UV and V-band luminosity function were in reasonable 
agreement with observations (Norberg et al. 2002; Kochanek et al. 
2001; Driver et al. 2012). 

Model origin The Galform-GP 14 model uses dark matter halo 
trees derived from the MS-W7 W-body simulation (Guo et al. 
2013), with a simulation box of 500 Mpc/fi side. 

Modifications to the supplied data In order to run Galform, we 
have remapped the given merger trees using D-Haloes (Jiang 
et al. 2014). This algorithm groups subhaloes in W-body cosmo¬ 
logical simulations avoiding transient structures and losses in mass. 
From the provided list of subhaloes, there is a percentage smaller 
than 15 per cent that D-Haloes classifies as independent haloes. 
This happens for haloes that when becoming subhaloes either re¬ 
tain more than 75 per cent their mass or that are located away from 
the main halo more than two half mass radii. 

Halo finder properties used The standard Galform models use 
Mbnd, but runs have been performed for all the five supplied mass 
definitions. 


A8 Munich - DLB07 (De Lucia & Blaizot) 

The variant of the Munich model described in De Lucia & Blaizot 
(2007b), with its generalization to the 3-yr Wilkinson Microwave 
Anisotropy Probe (WMAP3) cosmology discussed in Wang et al. 
(2008). The model includes prescriptions for gas cooling, star for¬ 
mation, stellar feedback, merger driven starburst, AGN feedback 
and chemical enrichment. The latter is based on an instantaneous 
recycling approximation. For more details on the physical models, 
we refer to Croton et al. (2006), De Lucia & Blaizot (2007b), and 
references therein. 

Cooling The rate of gas cooling is computed following the model 
originally proposed by White & Frenk (1991), and an implementa¬ 
tion similar to that adopted in Springel et al. (2001a). Full details 
can be found in De Lucia et al. (2010). 
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Star formation A Kennicutt-type prescription is adopted. Only gas 
above a critical surface density for star formation can be converted 
into stars. Details in Croton et al. (2006). 

Initial mass function A Chabrier (2003) IMF is used throughout. 

Metal treatment As detailed in De Lucia et al. (2004a), an instanta¬ 
neous recycling approximation is adopted. Metals are ejected (and 
instantaneously mixed) into the cold gas component after each star 
formation event. 

Supernova feedback and winds Supemovae explosions are as¬ 
sumed to reheat a (cold) gas mass that is proportional to the mass of 
stars formed. The amount of gas that leaves the dark matter halo in 
a wind is determined by computing whether the excess supernova 
energy is available to drive the flow after reheating the material to 
the halo virial temperature. Details can be found in Croton et al. 
(2006). 

Gas ejection & reincorporation Ejected gas is re-incorporated into 
the hot gas component on a time-scale that is related to the dynam¬ 
ical time-scale of the halo, as detailed in De Lucia et al. (2004a). 

Disk instability The instability criterion adopted is that of Efs- 
tathiou et al. (1982). When the instability condition is verified, we 
transfer enough stellar mass from the disk to the bulge so as to re¬ 
store stability. Eor details, see again Croton et al. (2006) and De 
Lucia et al. (2011). 

Starburst Galaxy mergers are accompained by starbursts mod¬ 
elled using the ‘collisional starburst’ prescription introduced by 
Somerville et al. (2001) with updated numerical parameters so as 
to fit the numerical results by Cox et al. (2004). 

AGNfeedback The model includes a distinction between a ‘quasar 
mode’ and a ‘radio mode’. AGN feedback is implemented as de¬ 
tailed in Croton et al. (2006). 

Merger treatment The model explicitly follows dark matter sub¬ 
structures. This allows us to follow properly the motion of the 
galaxies at their centres, until tidal truncation and stripping dis¬ 
rupt the subhalo at the resolution limit of the simulation. When this 
happens, a residual merger time is estimated using the current orbit 
and the classical dynamical friction formula of Binney & Tremaine 
(1987). Eor details, see De Lucia & Blaizot (2007b) and De Lucia 
etal. (2010). 

Substructures The model follows dark matter haloes after they are 
accreted onto a larger system, i.e. substructures are explicitly fol¬ 
lowed from the W-body simulation. 

Orphans An ‘orphan’ galaxy is created each time a substructure 
falls below the resolution limit of the simulation. The stellar mass 
of the galaxy is unaffected. The galaxy is assigned a residual 
merger time as detailed above. Its position and velocity are traced 
by following the most bound particle of the substructure at the last 
time it was present. As this information could not be reconstructed 
for this project positions and velocities of orphan galaxies are kept 
fixed to those of the substructures at the last time they were identi¬ 
fied. 


Calibration method The model has been calibrated ‘by-hand’ (note 
that in practice this means that a grid of model parameters was 
considered in order to evaluate the influence of each of them on 
model predictions). The main constrain is the K-band luminosity 
function at 2 = 0. 

Model origin The model is designed to work with merger trees 
from A-body simulations. 

Modifications to the supplied data No significant modification. 
Merger trees are reconstructed from original files provided, and us¬ 
ing each halo’s uniquely assigned descendant. 

Halo finder properties used The following quantities from the sup¬ 
plied halo catalogue are used: snapnum, positions and velocities of 
each halo, its mass M 200 C and spin parameter (this latter quantity 
is used to model the disk radius). 


A9 Munich - Lgalaxies (Henriques, Srisawat & Thomas) 

The ‘Munich’ model of galaxy formation is a semi-analytic scheme 
for simulating the evolution of the galaxy population as a whole 
and has been continually developed over the last quarter century 
(White & Erenk 1991; Kauffmann et al. 1993, 1999; Springel et al. 
2001b, 2005). Recent updates to the baryonic physics have resulted 
in a model that is capable of reasonably describing the observed 
population of galaxies in the local Universe. These include a de¬ 
tailed treatment of gas reheating, ejection and reincorporation by 
supernova (De Lucia et al. 2004b), updated in Guo et al. (2011), 
black hole growth during mergers (Kauffmann & Haehnelt 2000) 
and feedback from quiescent accretion (Croton et al. 2006), contin¬ 
uous environmental effects acting on satellite galaxies (Guo et al. 
2011), dust extinction from the ISM and molecular clouds (De Lu¬ 
cia & Blaizot 2007b). The model used in this paper, Henriques et al. 
(2013), includes all previous developments and aims at better rep¬ 
resenting the observed evolution of stellar mass across most of the 
age of the Universe. This is done by modifying the timescales for 
gas to be reincorporated after ejection from supernova in order to 
avoid an excessive build up of low mass galaxies at early times. As 
in Guo et al. (2013), a WMAP7 cosmology is adopted. 

In the latest major release of the Munich model, (Henriques 
et al. 2014), the modifications implemented in Henriques et al. 
(2013) were combined with a less efficient ram-pressure striping 
implementation in low mass groups and a lower threshold for star 
formation. These ensure that low mass galaxies are predominately 
star forming down to 2 = 0 and that the model can simultaneously 
match the evolution of the stellar mass function and the fraction of 
red galaxies. In addition, the AGN feedback and dust model were 
adjusted in order to better follow the properties of intermediate and 
high-z galaxies, respectively. Although the Henriques et al. (2014) 
model provides a significantly better representation of the observ¬ 
able universe, we are unable to run it for this project since it re¬ 
quires information on the trajectories of most-bound particles for 
haloes striped below resolution. These are used to follow the dy¬ 
namics of orphan galaxies and are crucially in order to track their 
properties. Without it, merger times and the disruption efficiency 
of satellites will be significantly different, changing the black-hole 
growth, star-burst efficiency and morphology evolution of all cen¬ 
tral galaxies. We therefore use Henriques et al. (2013) for which 
the combination of physics is less sensitive to these effects. 
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Cooling The cooling follows the implementation of White & Frenk 
(1991). Infalling diffuse gas is expected to shock-heat as it joins a 
halo. At early times and for low-mass haloes the post-shock cooling 
time is short and new material is assumed to settle onto the cold 
gas disk in a dynamical time. At later times and for higher mass 
haloes the shocked heated gas forms a quasi-static hot atmosphere 
from which it can gradually accrete to the centre via an element and 
temperature. 

Star formation Following Kauffmann (1996), star formation is as¬ 
sumed to be proportional to the mass in cold gas above a given 
threshold. 

Initial mass function A Chabrier (2003) IMF is assumed and 43% 
of the total mass of stars formed is assumed to be in massive, short 
lived stars, and immediately returned to the cold gas. 

Metal treatment Following De Lucia et al. (2004a), these stars are 
enrich the surrounding medium with a fixed yield of metals. Fu 
et al. (2013) and Yates et al. (2013) implemented, respectively, 
more detailed models of star formation and chemical enrichment 
in the Munich model. This will be incorporated in future releases. 

Supernova feedback and winds Supernova feedback is treated fol¬ 
lowing De Lucia et al. (2004a) and Guo et al. (2011). The fraction 
of the total available energy used in feedback is parametrised in a 
virial velocity dependent way. 

Gas ejection & reincorporation Part of the energy released during 
supernovae is used to reheat gas from the cold gas to the hot phase 
and the left-over used to eject material into an external reservoir. A 
new reincorporation model returns gas back into the hot phase at a 
rate inversely proportional to the halo mass (Henriques et al. 2013). 

Disk instability Disk instabilities are followed as in Guo et al. 
(2011) and transport material inwards to the bulge and they occur in 
galaxies where self-gravity of the disk dominates the gravitational 
effects of the bulge and halo. When the instability criteria is met, 
we transfer sufficient stellar mass from the inner parts of the disk 
to the bulge to make it marginally stable again. 

Starburst The stellar mass formed during a merger (see below) is 
modelled using the collisional starburst formulation of Somerville 
etal. (2001). 

AGN feedback Henriques et al. (2013) includes two black-hole re¬ 
lated processes: a “quasar” and a “radio” mode. Following Kauff¬ 
mann & Haehnelt (2000), black holes are taken to form and grow 
when cold gas is driven to the centre of merging systems. Black 
holes are also allowed to accrete gas from the hot gas atmospheres 
of their galaxies. This is assumed to generate jets and bubbles 
which produce radio mode feedback, suppressing cooling onto the 
galaxy and so eliminating the supply of cold gas and quenching star 
formation (Croton et al. 2006). 

Merger treatment When a satellite galaxy finally merges with the 
object at the centre of the main halo, the outcome is different for 
major and minor mergers. In a major merger, the disks of the two 
progenitors are destroyed and all their stars become part of of the 
bulge of the descendent, along with any stars formed during the 
merger. In a minor merger, the disk of the larger progenitor survives 


and accretes the cold gas component of the smaller galaxy, while 
its bulge accretes all the stars of the victim. Stars formed during the 
merger stay in the disk of the descendent. 

Substructures When the host halo of a galaxy enters the virial ra¬ 
dius of a larger system it becomes a satellite and its properties are 
strongly affected by a number of processes collectively called en¬ 
vironmental effects. Satellite galaxies do not receive primordial in¬ 
fall, they are stripped of their hot and ejected gas components due 
to ram-pressure and tidal striping, and they might merge with other 
galaxies (Guo et al. 2011). 

Orphans Once a satellite subhalo is disrupted, its central galaxy 
becomes an orphan and a merging clock is started. This will esti¬ 
mate how long the satellite will take to spiral into the central object 
due to dynamical friction. Since our implementation of tidal strip¬ 
ping of hot gas is directly connected to the stripping of dark matter, 
orphans have no hot or ejected gas. In addition, any cold gas re¬ 
heated by supernovae is added to the central halo and tidal forces 
might completely strip the stars and cold gas into the intracluster 
medium. 

In the default version of our model, the dynamical properties 
of orphan galaxies are given by those of the most bound particle 
identified at the time at which the halo falls below resolution. The 
vector offset is decayed due to dynamical friction. This has been 
shown to be crucial in order to correctly trace satellite distribu¬ 
tions and achieve convergence for simulations of different resolu¬ 
tion (Guo et al. 2011). Since the most bound particle information is 
not provided, the positions and velocities of orphans are frozen at 
the time of disruption and should be ignored. These do not have a 
significant impact on the physics in Henriques et al. (2013) which 
depend mostly on the independently calculated merger times but 
stop us from running the latest Henriques et al. (2014). 

Calibration method The best-fit model was chosen by fully sam¬ 
pling the allowed regions in parameter space using the MCMC 
methodology described in Henriques et al. (2009, 2013). The stel¬ 
lar mass, if-band and B-band Luminosity functions at z = 0,1, 2 
and 3 were used as observational constraints. 

Model origin Henriques et al. (2013) is built on A^-body merger 
trees following the method introduced by Springel et al. (2001a). 
The substructure are followed directly from the dark-matter simu¬ 
lation and the supplied data was used in unmodified form. The fol¬ 
lowing quantities from the supplied halo catalogue are used: snap- 
num, positions and velocities of each halo, its mass (M 2 ooc 200 ) 
and spin. 

Modifications to the supplied data The supplied data was used in 
unmodified form. We note that insufficient information was sup¬ 
plied to allow use of the latest version of L-Galaxies (Henriques 
et al. 2014) which better models the stripping of satellite galaxies 
within groups and clusters. 

Halo finder properties used The following quantities from the sup¬ 
plied halo catalogue are used: snapnum, positions and velocities of 
each halo, its mass (M 2 ooc 200 ) and spin. 

AlO Munich - Sage (Croton) 

The new Semi-Analytic Galaxy Evolution (SAGE) model is an up¬ 
dated version of that first presented in Croton et al. (2006). We only 
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highlight the significant changes here and point to the 2006 paper 
and Croton et al. (in prep.) for a full description of the rest of the 
model. 

Sage is publicly available through the Theoretical Astrophys- 
ical Observatory® (TAO; Bemyk et al. 2014), an online virtual labo¬ 
ratory that includes tools to add hundreds of magnitude filters to the 
galaxy output (with or without dust), construct custom light-cones, 
build images, and then download the mock data to your local ma¬ 
chine, all without any requisite programming knowledge. 

Cooling Cooling is handled as in Croton et al. (2006). An isother¬ 
mal sphere is assumed and a cooling rate estimated from a simple 
continuity argument. 

Star formation The SAGE model calculates the mass of cold gas in 
the disk that is above a critical surface density for star formation. 
New stars then form from this gas using a Kennicutt-type prescrip¬ 
tion. 

Initial mass function SAGE assumes a Chabrier (2003) IMF. 

Metal treatment SAGE follows the simplistic metal treatment in¬ 
troduced in De Lucia et al. (2004a). A yield of metals is produced 
from each star formation event and is recycled instantly back to the 
cold gas from very short-lived stars. 

Supernova feedback and winds Feedback from supernova in SAGE 
is a two step process. First, an assumed mass loading factor pushes 
cold gas out of the disk into the hot halo. Second, if enough energy 
from supernova has been added to the hot halo carried by this gas, 
some of the hot gas becomes unbound and is removed to an ejected 
reservoir. 

Gas ejection & reincorporation Gas can be ejected from the halo 
from supernova or quasar winds. Ejected gas can be reincorporated 
back into the hot halo at a rate in proportion to the dynamical time 
of the dark matter halo. 

Disk instability SAGE uses the Mo et al. (1998) approximation to 
determining when a disk becomes unstable. When so, enough exist¬ 
ing stars are transferred to the bulge to make the disk stable, along 
with any new stars as a result of a starburst. 

Starburst The SAGE model applies the collisional starburst model 
introduced in Somerville et al. (2001) to determine the mass of cold 
gas that becomes new stars during a merger. 

AGNfeedback SAGE uses the radio-mode AGN heating model in¬ 
troduced in Croton et al. (2006), and a new quasar-mode wind 
model introduced in Croton et al. (in prep.). 

Merger treatment Mergers are treated using the method described 
in Croton et al. (in prep.). Satellites are either merged with the cen¬ 
tral galaxy or added to the halo’s intra-cluster stars, depending on 
the subhalo survival time relative to an average expected based on 
its infall properties. 

Substructures Substructures are explicitly followed from the N- 
body simulation. 

® https://tao.asvo.org.au 


Orphans No orphans are used in SAGE. A decision as to the fate 
of a satellite galaxy has already been made and executed before its 
subhalo is lost below the resolution limit of the simulation. 

Calibration method SAGE is calibrated by hand using the 2 = 0 
stellar mass function, cold gas fraction, stellar metallicity-stellar 
mass relation, baryonic Tully-Fisher relation, and black hole-bulge 
mass relation. 

Model origin SAGE is an evolution of the Lgalaxies semi- 
analytic code which is based solely on A-body simulations. 

Modifications to the supplied data No modification were made to 
the supplied data. 

Halo finder properties used The primary mass used was M 200 C. 
Additional halo finder properties used to build galaxies are the peak 
value of the circular rotation curve, the position, and the spin pa¬ 
rameter. 


All HOD - Mice (Castander & Carretero) 

The MICE project^® is producing large simulations to help the de¬ 
sign and interpretation of large scale cosmological observational 
projects. In this paper, we use the galaxy mock generation code 
that has been developed within MICE. The galaxy mock code pop¬ 
ulates dark matter haloes with a hybrid Halo Occupation Distri¬ 
bution (HOD; ling et al. e.g., 1998; Scoccimarro et al. e.g., 2001; 
Berlind & Weinberg e.g., 2002) and SubHalo Abundance Match¬ 
ing (SHAM; Vale & Ostriker e.g., 2004; Tasitsiomi et al. e.g., 2004; 
Conroy et al. e.g., 2006) technique (Carretero et al. 2015; Castander 
et al. 2014; Crocce et al. 2013). Eollowing the HOD philosophy, 
we assume that haloes are populated by central and satellite galax¬ 
ies. We assign luminosities to central galaxies based on abundance 
matching taking into account the scatter between halo mass and lu¬ 
minosity. The HOD gives us the number of satellites in each halo. 
The satellite luminosities are drawn from the satellite luminosity 
function. We distribute satellites inside the haloes following a tri- 
axial ‘modified’ NFW, tweaked to match the observed clustering. 
The HOD parameters are also varied until we find an acceptable fit 
to the galaxy clustering as a function of luminosity. We assign ve¬ 
locities to the galaxies assuming a Gaussian velocity dispersion dis¬ 
tribution given by the halo mass (Bryan & Norman 1998). Lastly, 
we assign colours and SEDs with recipes that fit the clustering as a 
function of colour. We calibrate our method with local constraints 
given by the Sloan Digital Sky Survey (SDSS; York & et al. 2000) 
using the MICE Grand Challenge simulation (Fosalba et al. 2013; 
Crocce et al. 2013; Fosalba et al. 2013; Hoffmann et al. 2015) as 
starting halo catalogue. In particular we reproduce the galaxy lumi¬ 
nosity function (Blanton & et al. 2003; Blanton et al. 2005b), the 
colour-magnitude diagram (Blanton et al. 2005) and the SDSS clus¬ 
tering properties as a function of luminosity and color (Zehavi et al. 
2011). We extend our recipes to higher redshift applying evolution¬ 
ary corrections to the galaxy colours and then resampling from the 
cosmos catalogue (Ilbert & et al. 2009) galaxies with compatible 
luminosities and colours at the given redshift. 

Cooling — 

http://www.ice.cat/mice 
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Star formation We obtained the star formation rate from the dust- 
corrected UV flux of the galaxy SED. 

Initial mass function We assume a ‘diet’ Salpeter IMF (Bell & de 
Jong 2001). 

Metal treatment We compute the metallicity from the absolute 
magnitude using empirically determined relations. 

Supernova feedback and winds — 

Gas ejection & reincorporation — 

Disk instability — 

Starburst — 

AGN feedback — 

Merger treatment — 

Substructures We use a hybrid method to treat substructures. We 
can use the substructures provided by W-body, but if they are not 
available we can generate them analytically. 

Orphans We compute the expected number of satellites for each 
halo following an HOD prescription. If the halo contains fewer sub¬ 
haloes, we generate as many new satellites as the HOD predicts. 
We call these new satellites orphans in this context. We place them 
in the halo following a NEW profile with a concentration index 
expected for its halo mass. 

Calibration method The method has been calibrated to reproduce 
the galaxy luminosity function, the colour magnitude-diagram and 
the galaxy clustering as a function of luminosity and colour. The 
calibration is performed at low redshift and extrapolated at higher 
redshifts. The calibration has been done minimizing a where we 
have altered the input parameters manually. Note that our method 
has only been calibrated out to redshift z=1.5 and although we have 
computed quantities at higher redshifts in this paper, they are just 
an extrapolation of our recipies that we have not calibrated. So, take 
that in mind when trying to interpret the MICE results beyond that 
redshift. 

Model origin The model stems from W-body simulations, namely 
the MICE Grand Challenge W-body simulation. 

Modifications to the supplied data The supplied data has not been 
modified 

Halo finder properties used The base for our method are the halo 
masses and we adopted Mfof as the choice for the data presented 
here. We further use the following properties from the input cata¬ 
logue: haloid, hosthaloid, number of substructures, number of par¬ 
ticles, position, velocity, radius, peak value and position of the cir¬ 
cular rotation curve, velocity dispersion, and concentration. 

A12 HOD - Skibba (Skibba) 

The model used for this work is based on the halo model of galaxy 
clustering developed in Skibba et al. (2006) and Skibba & Sheth 


(2009). ‘Central’ and ‘satellite’ galaxy luminosities are modelled 
such that the luminosity function and luminosity dependence of 
clustering are the same as that observed. All galaxy properties and 
their occupation distributions are determined by halo mass, con¬ 
centration, and halo-centric position; therefore, all environmental 
correlations are entirely a consequence of the fact that massive 
haloes tend to populate dense regions (i.e., no ‘assembly bias’ is 
assumed). Satellites are assigned to subhaloes in the simulation, as¬ 
suming an abundance matching-like procedure that includes scatter 
between stellar mass, subhalo mass, and To include colours, 

we develop a prescription for how the colour-luminosity distribu¬ 
tion of central and satellite galaxies depends on halo mass, and this 
model is consistent with the observed colour mark correlation func¬ 
tion. I assume that the fraction of satellite galaxies which populate 
the red sequence increases with luminosity, and this fraction is in 
agreement with galaxy group catalogues (Skibba & Sheth 2009) 
and with the gradual quenching of satellites’ star formation Font 
et al. (2008). Stellar masses are based on luminosities and col¬ 
ors; since these are consistent with observations, it is not surpris¬ 
ing that the stellar mass function and clustering are consistent with 
observations as well. It is important to note that the model’s ob¬ 
servational constraints are robust only at M* > lO®/i“^M 0 and 
Afhaio > therefore, I have applied a mass threshold, 

below which this model should not be extrapolated or compared to 
other models (or to observations). 

We have made many updates and improvements to the model, 
which will be described in Skibba (in prep.). We are including 
colour and stellar mass gradients within haloes (Hansen et al. 2009; 
van den Bosch et al. 2008) and a dependence of the colour distri¬ 
bution on halo mass at fixed luminosity (More et al. 2011; Hearin 
& Watson 2013). In addition, the model includes a treatment of dy¬ 
namically unrelaxed systems, including some non-central bright¬ 
est halo galaxies, the fraction of which is constrained by SDSS 
and mock group catalogues (Skibba et al. 2011; Skibba & Maccio 
2011 ). 

The model provided here is a sort of HOD-SHAM hybrid, in 
which I populated subhaloes when possible and when an insuffi¬ 
cient number of resolved subhaloes were found, I distributed the 
satellites with my model’s standard prescription. I generously pop¬ 
ulated low-mass subhaloes in this model, and as a result the orphan 
fraction is relatively low. 

Cooling — 

Star formation The model currently only includes optical colours; 
a subsequent version will include a model of star formation rates. 
The model’s stellar masses apply a calibration from Zibetti et al. 
(2009). 

Initial mass function A Chabrier (2003) IMF is used throughout. 

Metal treatment Gas-phase and stellar metallicities are based on 
scaling relations of Tremonti et al. (2004) and Gallazzi et al. (2005). 

Supernova feedback and winds — 

Gas ejection & reincorporation — 

Disk instability Disk instabilities are not included in the model, but 
spiral and elliptical morphologies are included based on clustering 
and other constraints from Skibba & Sheth (2009). 
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Starburst — 

AGN feedback AGN feedback is not modelled here, but a simple 
black hole mass scaling relation basted on Tundo et al. (2007) is 
applied. 

Merger treatment — 

Substructures Substructure properties are taken from the simula¬ 
tion at a given snapshot, but subsequent subhalo evolution is not 
modelled. 

Orphans When a sufficient number of substructures are not re¬ 
solved to match them to satellites, haloes are populated with the 
remaining satellites (‘orphans’) in order to reproduce the model’s 
occupation distributions. Orphans are spatially distributed with a 
Navarro et al. (1997, NFW) profile with the Maccio et al. (2008) 
mass-concentration profile is assumed, while accounting for the 
fact that galaxies and subhaloes are typically less concentrated than 
dark matter Munari et al. (2013). 

Calibration method The model has been designed to reproduce 
real- and redshift-space luminosity and colour-dependent galaxy 
clustering (Zehavi et al. 2011) and mark clustering statistics 
(Skibba et al. 2006). The model is also constrained by the SDSS 
luminosity function (Blanton &. et al. 2003; Yang et al. 2009) and 
colour-luminosity distribution (Skibba & Sheth 2009). The model 
is consistent with the Moustakas et al. (2013) stellar mass function, 
which is not used as a constraint. The model is consistent with the 
Moustakas et al. (2013) stellar mass function, which is not used as 
a constraint though. 

Model origin The model originates from an analytic halo-model 
formalism (Cooray & Sheth 2002). 

Modifications to the supplied data No modification were made to 
the supplied data. 

Halo finder properties used The model uses the 3-D positions, 3-D 
velocities, halo mass and radius (200c and bound), Vmax, velocity 
dispersion (200c), and the substructure abundances and properties. 

Because the observational constraints are less robust at low 
masses and luminosities, we have applied a halo mass threshold 
near The SKIBBA model is complete only above this 

mass. 


APPENDIX B: HALO-MASS DEFINITION 

While we have already discussed the influence of the applied mass 
definition on the model-to-model variation seen in Section 5 we like 
to extend this here a bit more by directly comparing two mass def¬ 
initions on a model-by-model basis. Namely, for those models that 
provided both an M 200 C galaxy catalogue and a catalogue based 
upon their own mass definition (different from M 200 C, of course) 
we show in Fig. B1 the stellar mass functions: solid lines are for the 
model mass definition whereas dotted lines are for M 200 C. Some 
models appear to be sensitive to the choice of the definition for 
halo mass, but the overall level of scatter across various models 
seems to stay similar. We append some more plots that focus on 



Figure Bl. Stellar mass function at redshift z = 0 for models that also 
provided data for M 200 C as the mass definition. To be compai'ed against 
Fig. 2. 



Figure B2. Stellar-to-halo mass ratio for all (non-oiphan) galaxies at red- 
shift 2 = 0. To be compared against Fig. 10. 

the redshift evolution as this is where the mass definition leaves its 
largest imprint. The plots are accompanied by Table B1 where we 
list the number of galaxies in various populations. This table should 
be compared against Table 3. 


APPENDIX C: INITIAL STELLAR MASS FUNCTION 

As the change in stellar mass will also influence the number of 
galaxies above our usual threshold we also list the 

(change in) numbers in Table Cl, only showing the affected mod¬ 
els. This table should be compared against Table 3 again. 


APPENDIX D: UN-NORMALIZED REDSHIFT 
EVOLUTION 

In Section 6 we discussed the redshift evolution of both the number 
(density) of galaxies and the star formation rate (density), normal¬ 
izing the respective curves to their redshift 2 = 0 values. The nor¬ 
malizations have been provided in Table 2 & Table Cl and hence 
we separated ‘trends’ (as shown in the figures) from ‘absolute’ 
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Figure B3. Star formation rate density as a function of redshift for models 
that also provided data for M200C as the mass definition. To be compai'ed 
against Fig. 7. 



Figure B4. The number density of all galaxies with stellar mass M* > 
as a function of redshift for models that also provided data for 
M 20 OC as the mass definition. To be compared against Fig. 11. 

Table Bl. Number of galaxies at redshift 2 = 0 with a stellar mass in excess 
of M* > for models when applying M20OC as the mass def¬ 

inition (—200c extension) and when using their favourite mass definition. 
Note again that for MORGANA all satellites have been tagged as ‘orphan’. 


Table Cl. Number of galaxies at redshift 2 = 0 with a stellar mass in 
excess of M* > lO^/i“^M0 for those models that have been affected by 
the transformation to a Chabrier IMF. To be compared against Table 3. 


code name -^gal -^central -^non — orphan -^orph an 

GalICS- 2.0 tuned parameters to observations w/ Chabrier IMF 


Sag 

14025 

9565 

11758 

2267 

Durham flavours: 

Galform-gp14 

9842 

5680 

6770 

3072 

Galform-kb06 

10467 

6060 

7225 

3242 

Galform-kf08 

13340 

7054 

8371 

4969 

HOD models: 

Mice 

9510 

5638 

8067 

1443 



Figure Dl. Star formation rate density as a function of redshift. 


differences (as listed in the tables). Here we now provide the un¬ 
normalized plots directly showing the different evolutions for both 
these quantities. 

This paper has been typeset from a TgX/ HTeX file prepared by the 
author. 


code name 

gal 

-^central 

N^=o 

non —orphan 

Nz=o 

orphan 

Galacticus 

14255 

7825 

10019 

4236 

Galacticus— 200c 

16123 

9026 

11393 

4730 

Galform-gp14 

8824 

5097 

6098 

2726 

Galform-gp14-200c 

9320 

5595 

6666 

2654 

Morgana 

10008 

6186 

6186 

3822 

Morgana— 200c 

7316 

4925 

4925 

7316 

Sag 

19516 

13571 

16256 

3260 

Sag-200c 

16505 

11332 

13773 

2732 
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Figure D2. The number of all galaxies with stellar mass M* > 
as a function of redshift. 
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